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ABSTRACT 

Optimal  age  replacement  policies  are  designed  to  cut  down  system  failures  and 

minimize  maintenance  cost.  By  scheduling  planned  replacements,  a  system  is  replaced 
at  age  <f>*  or  at  failure,  whichever  comes  first,  and  the  cost  of  replacement  before  failure 
(planned)  is  less  than  the  cost  after  failure  (unplanned).  In  this  thesis,  the  distribution  of 
lifetimes  is  a  known,  increasing  failure  rate  phase  type  distribution.  To  find  the  optimal 
age  of  replacement,  the  parameters  of  the  underlying  phase  type  distribution  must  be 
estimated. 

An  optimal  age  sequential  estimation  procedure  is  developed.  In  particular,  the 
phase  type  distributions  parameters  are  estimated  using  a  matching  moments  nonlinear 
programming  approach.  Since  there  are  many  parameters  associated  with  phase  type 
distributions  and  the  distributions  include  matrix  exponential  terms,  the  parameters  are 
in  general  difficult  to  estimate.  A  specific  case  where  the  phase  type  distribution  has 
initial  probability  vector  a =(1,0,0)  is  studied  for  different  sample  sizes  and  compared 
with  a  similar  nonparametric  procedure. 
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I.   INTRODUCTION 

Normally  a  system  or  component  is  replaced  when  it  fails  and  Cu  the  cost  of 
replacement  before  failure  (planned  replacement),  is  often  less  than  C2,  the  cost  after 
failure  (unplanned  replacement).  The  problem  of  determining  when  to  repair  and  when 
to  replace  failing  systems  is  the  concern  of  management  resources.  Inefficient 
management  due  to  the  use  of  nonoptimal  maintenance  policies  can  lead  to  significant 
system  maintenance  cost.  In  general,  optimal  maintenance  policies  are  designed  to  cut 
down  the  number  of  system  failures  and  minimize  maintenance  costs.  In  this  thesis,  we 
study  the  problem  of  estimating  a  type  of  optimal  maintenance  policy  (the  optimal  age 
replacement  policy)  when  the  underlying  system  life  distribution  is  phase-type.  In 
particular,  we  consider  an  adaptive  estimation  scheme  where  the  estimated  optimal 
replacement  age  is  updated  each  time  the  system  is  replaced. 

Maintenance  is  defined  to  be  all  activities  taken  to  keep  the  system  in  serviceable 
condition  or  to  bring  it  back  to  serviceability.  There  are  two  types  of  maintenance, 
corrective  and  preventive  maintenance  [Ref.  5:  pp.  1-8].  Corrective  maintenance  is  used 
after  a  failure.  This  does  not  necessarily  mean  that  such  action  has  not  been  foreseen. 
Preventive  maintenance  aims  to  reduce  the  probability  of  failure  and  includes: 

-Planned  (scheduled)  maintenance  in  which  specified  components  are  replaced 
(usually  at  regular  intervals).  The  maintenance  time  is  usually  based  on  component 
lifetime  distributions. 


-Unplanned  (condition-based)  maintenance  in  which  the  decision  to  replace  or  not 
to  replace  is  made  according  to  the  outcome  of  a  diagnostic  study. 

The  policies  that  are  designed  to  reduce  the  number  or  the  probability  of  system 
failure  and  maintenance  costs  are  called  maintenance  policies  [Ref.  21:  pp.  1-3].  We 
concentrate  on  age  replacement  policies  where  a  system  is  replaced  at  age  T  or  failure, 
whichever  comes  first.  For  this  problem  to  make  sense,  Q  must  be  less  than  C2  and  the 
system  must  age  with  time  i.e.  the  failure  rate  of  the  system  must  be  increasing  with  age. 
It  is  not  wise  to  replace  the  equipment  before  failure  when  the  system  has  a  constant  or 
decreasing  failure  rate  such  as  when  the  system  failures  occur  according  to  an  exponential 
distribution.  In  this  case,  replacement  will  not  reduce  the  probability  of  system  failure 
occurring  in  the  next  instance  of  time.  When  using  an  age  replacement  policy  the 
question  always  asked  is,  "At  what  age  should  the  equipment  be  replaced".  If  the 
replacement  occurs  too  frequently  the  maintenance  costs  will  be  high.  If  replacement  is 
too  infrequent,  the  system  will  fail  more  often  than  necessary  and  again,  the  maintenance 
cost  will  be  too  high.  Thus  it  is  desirable  to  find  an  "optimal"  replacement  age.  Here, 
the  optimal  age  is  defined  as  the  age  which  yields  the  minimum  long  run  expected 
maintenance  costs  per  unit  time. 

We  will  assume  that  the  maintenance  action  returns  the  system  to  the  good-as-new 
condition,  thus  the  same  services  are  provided  as  before  replacement.  To  accomplish  this 
we  assume  that  the  systems  used  for  replacement  have  lifetimes  that  are  independent  and 
identically  distributed  according  to  a  phase-type  distribution.  The  class  of  phase  type 
distributions  is  large  and  includes  Exponential  and  Gamma  distributions  along  with 


convolutions  and  mixtures  of  these  distributions.  This  fact,  along  with  the  fact  that  phase 
type  distributions  have  a  physical  interpretation  make  them  particularly  well  suited  for 
modeling  system  lifetimes. 

The  optimal  replacement  age  depends  on  the  underlying  phase  type  distribution. 
This  is  in  general  unknown  and  must  be  estimated.  Estimation  for  phase  type 
distributions  based  on  iid  lifetimes  has  been  studied  by  Neuts  and  Meier  (1980). 
Estimation  of  the  optimal  age  of  replacement  for  phase  type  distributions  is  new. 
However  both  parametric  and  nonparametric  estimation  based  on  iid  observations  have 
been  considered  by  Arunkumar  (1972),  Ingram  and  Schaeffer  (1976),  Bergman  (1977), 
and  Barlow  (1978).  Here  we  use  a  sequential  approach  for  estimating  the  optimal 
replacement  age.  At  each  replacement  the  estimate  is  updated  and  the  new  system  is 
subject  to  an  age  replacement  policy  based  on  the  optimal  age  estimated  so  far.  This  type 
of  sequential  approach  has  been  studied  parametrically  by  Oclay  (1990)  and 
nonparametricly  by  Bather  (1977),  Frees  and  Ruppert  (1985),  Aras  and  Whitaker  (1991), 
Wu  (1990)  and  in  a  decision  theoretic  framework  by  Glazebrook,  Bailey  and  Whitaker 
(1991). 

Phase  type  distributions  are  described  in  Chapter  II .  The  sequential  estimation 
procedure  to  estimate  the  optimal  age  of  replacement  is  given  in  Chapter  HI,  and  in 
Chapter  IV  we  present  results  comparing  the  sequential  estimation  procedure  assuming 
an  underlying  phase  type  distribution  with  a  nonparametric  procedure.  Conclusion  and 
recommendation  are  given  in  Chapter  V. 


n.   PHASE  TYPE  DISTRIBUTION 

A.      GENERAL 

A  phase  type  distribution  is  defined  as  the  distribution  of  the  time  until  absorption 
in  a  finite  state  Markov  process.  To  determine  the  distribution  of  the  absorption  time, 
consider  anm+1  state,  continuous  time  Markov  chain  whose  infinitesimal  generator  Q 
has  the  form 


Q    = 


T  T' 

0    0 


where  T  is  a  nonsingular  m  x  m  matrix  with  negative  diagonal  elements  and  nonnegative 
off-diagonal  elements,  T°  is  vector  of  length  m  with  nonnegative  elements  and 
0 = (0,0,..., 0).  Moreover, 

Te  +  T°  =  0' 
where  e  =  (1,1,.. .,1)'. 

By  the  construction  of  Q  the  state  m+1  is  absorbing  and  all  other  states  are 
transient.  The  necessary  and  sufficient  condition  for  this  is  that  the  inverse  T1  exists 
[Ref.  15:  pp.  446].  Let  the  initial  probability  vector  of  the  Markov  chain  be  given  by  ( 
«>  «m+i )  with  «e  +  o^+x  =  1  and  a  being  the  m  dimensional  initial  probability  vector 
of  transient  states  such  that  0  <  ae  ^  1 .  Let  X  be  the  time  until  absorption  into  the 
(m+l)rt  state.  The  probability  distribution  F(x)  of  the  time  until  absorption  in  the  state 
m+1  corresponding  to  the  initial  probability  vector  (a,  c^+i )  is  given  by 


F(x)  =  1  -  aexp(Tx)e  ,      for  x  £  0  ,  (2.1) 

where  exp(A)  is  the  matrix  exponential  of  the  matrix  A  defined  in  Ref.  8. 

The  probability  distribution  F(x)  on  [0,  oo)  is  said  to  be  of  phase  type  (PH- 
distribution)  and  the  pair  (a,  T)  is  called  a  representation  of  F(x).  If  c^+i  >  0  the 
distribution  F(x)  has  a  jump  of  height  c^+i  at  x  =  0. 
On  (0,  oo ),  the  probability  mass  is  distributed  according  to  a  density  given  by 

f(x)  =  oexp(Tx)T°  ,  for  x  >  0  ,  (2.2) 

The  noncentral  moments  /x;  =  E[X']  of  F(x)  are  all  finite  and  given  by 

K  =  (-lyiKoTe)  ,   for  i  2:  1  [Ref.  15:  pp.  446]  .  (2.3) 

B.      COST  FUNCTION 

Let  Xlt  X2, . . .  be  a  sequence  of  independent  and  identically  distributed  (iid)  positive 
random  variables  with  distribution  function  F.  The  sequence  Xlf  X2,...  represents  the 
sequence  of  system  lifetimes  that  would  be  observable  if  the  system  were  replaced  at 
failure.  Let  Cl5  Q  (C2  >  Q)  and  <i>  be  the  respective  cost  of  planned  replacement; 
unplanned  replacement  and  the  replacement  age.  The  long  run  expected  cost  per  unit  time 
can  be  verified  [Ref.  1:  pp.  87]  to  be 

_     c2f(4>)  +  qF(») 

C(<P}      "     J^Z (2.4) 

\    F(u)du 
Jo 

where  F(u)  =  1  -  F(u)  is  the  survival  function.  A  sufficient  condition  that  guarantees  a 
unique  and  finite  <f>'  that  minimizes  C(<£),  is  that  F  have  strictly  increasing  failure  rate. 


When  the  distribution  F  is  phase  type  with  representation  (a,  T),  the  long  run  expected 
cost  function  is 


C(<J>)      = 


C2  [l-aexp(2V|>)  e]    +  Cx[aexp(3\|>)e] 


f  aexp(Tu) edu 
Jo 


(2.5) 


C2  -  ttexp(a^>)e[C2-C1] 
oT_1[exp(!ntJ))    -  I]  e 


(2.6) 


where  I  is  an  identity  matrix. 

The  phase  type  distribution  does  not  necessarily  have  increasing  failure  rate,  so  the 
cost  function  C(<j>)  does  not  necessarily  have  a  unique  finite  minimum.  In  the  m =3  case 
even  when  the  initial  probability  of  the  transient  states  is  fixed  the  minimum  can  be 
unique  and  finite  or  can  occur  at  infinity  as  shown  in  Figure  1,  where  a  =  (1,  0,  0)  and 
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In  Figure  2,  the  same  phenomena  is  shown  when  F  has  the  same  nonsingular  3x3  matrix 
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but  the  initial  probabilities  are  different. 
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Figure  1  The  long  run  expected  cost  per  unit  time  curves  of  PH  distribution  with 
the  same  a  and  different  T  =  T,  or  T2 
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Figure  2  The  long  run  expected  cost  per  unit  time  curves  of  PH  distribution  with 
the  same  T  and  different  a 
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m.    THE  SEQUENTIAL  ESTIMATION  PROCEDURE 

In  theory,  any  distribution  of  a  nonnegative  random  variable  can  be  approximated 
arbitrarily  closely  by  a  PH  distribution  [Ref.  11:  pp.  1-2].  This  means  that  the  PH 
distributions  are  dense  in  &~  where  ^"is  the  set  of  distributions  with  support  on  [0,oo). 
The  paper  of  Johnson  and  Taaffe  (1988)  shows  that  the  finite  mixtures  of  Erlang 
distributions  and  PH  distributions  with  a  Coxian  representation  are  both  dense  in  &.  Since 
both  families  are  subsets  of  the  family  of  PH  distributions,  this  implies  that  PH 
distribution  are  dense  in  ^[Ref.  10:  pp.  1-8].  Thus,  the  PH  distribution  can  be  used  to 
approximate  any  unknown  distribution  of  lifetimes  (lifetime  is  always  greater  than  or 
equal  0).  Another  feature  that  makes  PH  distributions  a  desirable  choice  to  model  system 
lifetimes  is  that  they  may  have  a  physical  interpretation.  The  absorbing  state  represents 
system  failure  and  the  transient  states  may  represent  different  levels  of  a  system's  ability 
to  function. 

A.      PARAMETER  ESTIMATION 

There  are  many  ways  to  estimate  the  parameters  of  a  distribution  including  the 
method  of  maximum  likelihood  (MLE),  and  the  method  of  moments.  The  PH 
distributions  have  special  properties  which  make  estimation  difficult.  First,  the  number 
of  parameters  is  not  fixed,  it  varies  with  the  number  of  transient  states,  e.g.  a  PH 
distribution  with  2  transient  states  has  6  parameters,  a  PH  distribution  with  3  transient 
states  has  12  parameters  etc.  Also  the  probability  distribution  and  density  function  include 


the  exponential  of  a  matrix.  This  makes  it  hard  to  work  with  the  likelihood  function.  In 
this  thesis,  the  moment  matching  method  is  used  to  estimate  the  parameters  of  a  phase 
type  distribution  but  the  MLE  method  and  its  associated  problems  are  also  discussed  . 
In  addition,  the  number  of  transient  states,  m,  is  assumed  to  be  known. 

1.       The  Maximum  Likelihood  Method 

Let  Xu  X2,...  represent  the  sequence  of  system  lifetimes,  where  Xls  X2)... 
are  iid  PH  distributions  with  representation  (a,  T).  After  N  observations  the  data 
available  for  estimation  will  be  (Zl6l)t(Z2,dd,...,(ZK,bN)  where  Z^min  (&,$*&),  X{  is 
the  ^lifetime  and  $*iA  is  the  optimal  age  replacement  estimated  after  i-1  observations  and 
5j=l  if  Xi<$*iA  and  8—0  otherwise.  Assuming  system  lifetimes  have  a  PH  distribution 
the  likelihood  for  the  replacement  ages  Zj,...,^  and  types  of  replacement  8U  62,...,6n 
is 


L(a,T)    =  IJ  [aexpCrz^r'l^taexpdTZ^e]1"^ 
i=i 


IT 

=  IJ  [-oexp(rzi)re]*i[aexp(rzi)e] 
1=1 


It  is  too  difficult  to  maximize  this  likelihood  directly  by  differentiating  L(a,T)  and 
solving  for  the  MLE's.  There  are  several  reasons  for  this.  The  likelihood  and  its 
derivatives  include  the  exponential  of  a  matrix  form  that  cannot  be  written  in  closed  form 
and  the  parameters  are  subject  to  numerous  constraints.  In  addition,  there  are  many 
likelihood  equations  due  to  the  number  of  parameters  (3  transient  states  will  have  up  to 
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12  equations).  An  alternate  approach  is  to  find  approximate  MLE's  using  nonlinear 
programming  algorithms  i.e.: 

m 
MAXL(a,T)    =  JJ  [-oexp  (TZ±)  Te]  &i  [aexp  (TZj)  e] 1_8i 
i=i 

S.T.  a         ^0 

ae      =     1 

^      <     0  fori  =  1,2,.  ..,m 

^      £     0  for  i  *  j 

Ei^tjj      <    -t,  fori  =  1,2,. ..,m 

det(T)     *    0 

where  m  =  the  number  of  transient  states  of  the  PH  distribution.  Even  when  m  is 

assumed  to  be  known,  there  is  still  the  problem  of  approximating  the  exponential  of  a 

matrix  and  the  optimization  software  to  take  care  of  this  problem  is  not  available.  Thus, 

it  is  very  difficult  to  use  this  approach  to  estimate  the  PH  distribution  parameters. 

2.       The  Moment  Matching  Method. 

The  PH  distribution  is  a  complicated  distribution,  there  are  many  parameters 
and  there  are  difficulties  with  computing  the  exponential  of  a  matrix.  One  property,  the 
existence  of  T1,  implies  that  all  noncentral  moments  of  a  PH  distribution  are  finite.  The 
k*  moment  is  given  by 

Aik   =    (-l^la^e.  (3.1) 

Moment  matching  is  a  common  method  for  estimating.  Using  moment  matching  by- 
passes the  problem  of  evaluating  the  exponential  of  a  matrix  [Ref.  12:  pp.3]. 
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Let  m  be  the  number  of  transient  states  of  the  PH  distribution.  There  are  m(m+l) 
parameters  that  need  to  be  estimated,  thus  the  m(m+l)  equations  from  the  moments  that 
need  to  be  solved  are 

£,     =    -aT^e 

p,2     =   2!aT1e 


/Wd    =    (-ir(ffl+V+m)!«A 
where  /xk  is  the  k*  sample  moment. 

This  method  is  still  inappropriate  due  to  the  large  number  of  equations  and 
the  fact  that  when  the  dimension  of  matrix  T  is  big  and  the  moments  are  of  high  order, 
substantial  amount  of  error  is  introduced  when  trying  to  solve  these  equations.  Johnson 
and  Taaffe  have  shown  that  the  moment  matching  method  and  nonlinear  programming 
approaches  can  be  used  to  estimate  the  parameters  of  a  PH  distribution.  They  match  only 
the  second  and  third  standardized  moments  [Ref.  11:  pp.  3-11]. 

Let  fij  be  the  i*  central  moment.  The  second  standardized  moment,  the  coefficient 
of  variation  (C),  defined  as 

c    _     standard  deviation 
mean 

can  be  written  as 
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C    =       -J-2- .  (3.2) 


The  third  standardized  moment,  the  coefficient  of  skewness,  (7)  is  defined  as 


(3.3) 


(TT2)3/2 

And  the  t*  standardized  moment  is  fij(fx2)t'2  for  t  =  3,4,...  . 

The  reason  for  matching  standardized  moments  is  that  the  standardized  moments  do  not 
change  with  scale  changes.  So,  the  shape  of  a  PH  distribution  with  representation  (a,T) 
will  not  change  if  it  is  multiplied  by  a  nonnegative  number. 

The  nonlinear  programming  formulation  for  matching  C  and  7  to  a 
continuous  PH  distribution  is  given  by 

MIN.  w,(C-C(0))2  +  w2  (7-7(0))2 

S.T. 

a    >    0 
ae  =    1 

fa  <    0  fori  =  l,2,...,m 

tg  £    0  for  i  *  j 

Ei*jtij  ^  -ta  fori  =  1,2,. ..,m 

det(T)  *  0 
where 

C(0)  =  the  second  standardized  moment  of  the  current  solution. 
7(0)  =  the  third  standardized  moment  of  the  current  solution. 
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m  =  the  number  of  transient  states  of  the  PH  distribution. 
tjj  =  the  elements  of  row  i  and  column  j  of  nonsingular  matrix  T. 
wl5w2  are  positive  weights  chosen  to  guide  the  search  and  wt  ^  3w2. 
The  moments  are  matched  when  the  objective  function  value  is  zero.  The  solutions  which 
match  the  target  moments  are  called  "moment-matching  solutions"  [Ref.  11:  pp.  5]. 

Even  with  fixed  dimension,  this  procedure  may  have  many  solutions  and  these 
solutions  are  unpredictable.  There  are  some  practical  points  that  need  to  be  made 
associated  with  the  properties  of  PH  distribution  and  the  problem  of  using  nonlinear 
programming  to  get  the  moment  matching  solutions: 

-  Different  combinations  of  initial  solutions  and  algorithm  parameters  may  lead  to 
different  moment-matching  solutions.  Also,  some  combinations  of  initial  solution  and 
algorithm  parameters  may  not  lead  to  a  moment-matching  solution  and  the  nonlinear 
programming  algorithm  may  not  converge. 

-  For  a  given  dimension  we  do  not  know  how  many  solutions  exist  or  how  to  guide 
the  search  toward  a  preferable  solution. 

-  The  moment-matching  solutions  do  not  guarantee  that  the  estimated  cost  function 
will  have  a  finite  minimum.  When  choosing  among  several  solutions,  a  solution  which 
gave  a  cost  function  with  a  finite  minimum  was  always  chosen. 

-  The  values  of  W!  and  w2  can  be  modified  to  guide  the  search;  generally  w,  ^  1, 
w2^l  and  W!  ^  3w2.  The  magnitude  of  w,  and  w2  can  be  increased  to  get  more  or 
sufficient  accuracy  [Ref.  11:  pp.  6]. 
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-  If  the  feasible  solutions  are  known,  parameters  can  be  bounded  to  guide  the 
search.  However,  some  bounds  on  the  parameters  may  not  lead  to  a  moment-matching 
solution  or  to  a  solution  whose  cost  function  does  not  have  a  unique  finite  minimum. 

-  User  interaction  is  often  necessary  in  obtaining  the  appropriate  or  preferable 
solutions. 

B.      THE  SEQUENTIAL  ESTIMATION  PROCEDURE 

Let  XUX2, ,XN  be  the  sequence  of  lifetimes  of  the  system  and  {$*N}  be  the 

sequence  of  estimators  of  <f>*  where  $*N  is  the  estimator  after  N  replacements.  Then  after 
N  observations,  we  have  the  right  censored  data  (Z1,51),(Z2,52),....,(ZN,6N). 
where         Z{  =  min  (Xj,$*M) 

Si  =  1  if  X;  ^  $Vi  otherwise  8{  =  0  . 
Because  the  data  are  right  censored  the  usual  method  of  moments  approach  for  estimation 
needs  to  be  modified.  Rather  than  use  sample  moments  calculated  from  an  empirical 
distribution,  we  use  moments  from  a  nonparametric  estimator  of  F  based  on  the  censored 
data.  From  the  right  censored  data  after  N  observations,  the  procedure  to  compute  the 
estimators  {$*{}  is  developed  as  follows. 

1.  Use  the  right  censored  data  sequence  (Z^^),^,^),  ...,("ZnM  to  estimate 
moments  with  F  =  l-F  by  the  product-limit  estimator: 
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*W      =  II         (  ^~  )8(i>  (3.4) 

where  Z(1),  Z^,...,  Z^  are  the  order  statistics  of  Z,,  Z2,...,  ZN  and  5(1),  6(2),...,  5^  are 
ordered  according  to  the  ordering  of  Z(1),  ZQ),...,  Z^.  Then  calculate  probabilities  P;, 
associated  with  Z(i),  by 

P(  =  F(Z(i)")  -  F^)  .  (3.5) 

Note,  since  F  has  discontinuities  only  at  Z&  where  5(i)  =  1 ,  that  P(  =  0  when  5(0  =  0. 
Then  estimates  of  the  k*  moment  are  given  by 


£[xk]     =    J2pizku)  (3'6) 

The  estimate  of  the  2nd  standardized  moment  is 

£      _      y/£[X2]  ~(E[X]  )2  (3    7) 

E[X] 

and  the  estimate  of  the  3rd  standardized  moment  is 

f       =     g[X3]  -3E[X]E[X2]  +  2(E[X])3        {3    Q) 
(y/E[X2)  -(EiX])2)* 

2.  Estimate  parameters  of  PH  distribution  (a,T),  by  matching  the  second  and  third 
standardized  moments  using  the  nonlinear  programming  approach. 

When  executing  the  nonlinear  programming  algorithm  using  the  model  described 
in  the  previous  section,  the  det(T)  ^  0  constraint  is  replaced  by  a  constraint  on  the 
expected  lifetime.  The  reason  for  doing  this  is  that  we  could  not  formulate  an  algebraic 
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constraint  equivalent  to  the  constraint  det(T)  *  0  .  The  problem  can  be  solved  by  a 
constraint  on  the  expected  lifetime,  i.e.  by  taking 

oT"1©    =    £[X]  (3.9) 

where  E[X]  is  calculated  in  step  1 .  This  constraint  will  make  the  search  for  preferable 
solutions  easier.  The  initial  solution  is  chosen  by  the  user.  Some  initial  solutions  may  or 
may  not  lead  to  a  moment  matching  solution.  Typically,  the  initial  vector  a  is 
(l/m,l/m,...,l/m)or(l,0,...,0)and  an  initial  matrix  T  consists  of  tM  =  -1  and  tjj  =  1/m 
for  i  ?*  j,  i  =  l,2,...,m.  In  this  thesis  the  initial  vector  a  is  (1,0,..., 0)  and  an  initial 
matrix  T  is  taken  to  have  t;i  =  -1  and  t^  =  1/m  [Ref.  11:  pp.  9]  . 

3.  Using  the  parameters  estimated  in  step  2,  the  estimated  cost  function  is 

C(*)     =     C2-aexp(^)e[C2-C1]  (3   ±Q) 

oT^texpffty)    -  l]e 

The  new  estimate  of  the  optimal  age  of  replacement  $'N  is  taken  to  be  the  <f>  that 
minimizes  (3.10)  by  enumerative  search. 

4.  Compare  the  $*N  with  XN+1  then  repeat  Step  1. 

The  initial  solutions  for  matrix  T  and  vector  a  in  step  2  are  taken  to  be  the  previously 
estimated  values  of  T  and  a.  In  case  the  moment  matching  solution  cannot  be  met,  the 
initial  solution  in  step  2  will  be  taken  to  be  the  original  initial  solution  and  the  previous 
optimal  age  replacement  is  used  for  step  3. 

The  replacement  cost  for  i*  system  is  C2  if  X;  <  0"M;  otherwise  the  replacement 
cost  is  Ct.  The  actual  total  replacement  cost  for  the  first  N  systems  that  are  observed  is 
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N 


Cn    =    E  [C2x6i  +  qx  (1-6^]  (3.11) 

i=l 

and  the  total  operating  time  for  the  N  systems  is 

tN    =    f^MiniXt.fr^)     -    fjz,    .  (3-12) 

i=l  i=l 

So,  the  actual  average  cost  (AAC)  after  N  replacement  is  computed  by 


AACN    =     -^  (3.13) 

In  this  thesis,  the  moment  matching  nonlinear  programming  approach  for  estimating 
parameters  of  PH  distributions  uses  the  GAMS  program  as  shown  in  Appendix  A  (for 
a  PH  distribution  with  3  transient  states).  The  Fortran  programs  are  written  to  estimate 
the  optimal  age  of  replacement,  second  and  third  standardized  moments  and  for  data 
preparation  as  shown  in  Appendix  B. 
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IV.COMPARISON  WITH  NONPARAMETRIC  PROCEDURE 

A.      NONPARAMETRIC  PROCEDURE 

An  alternate  approach  for  estimating  <t>*  is  to  estimate  nonparametrically,  as  in  Aras 
and  Whitaker  (1991).  The  procedure  is  much  simpler  computationally  than  the  parametric 
procedure  described  in  the  previous  section  for  PH  distributions.  In  the  nonparametric 
procedure  after  N  replacements  the  product  limit  estimator  F  of  F  given  in  (3.4)  is  used 
to  estimate  the  cost  function  C(<f>)  as 

r  <*>)    -    c2f(<|>)  +  q#(j>) 

cnW>      "     ^ (4.1) 

[  F(u)du 
Jo 

where  F  =  1  -  F  is  the  estimator  of  the  cdf  F,  the  estimator  $*N  of  <f>'  is  then  found  by 

minimizing  CN(<£).  In  general  one  would  expect  parametric  procedures  to  do  much  better 

than  nonparametric  procedures.  However,  with  the  numerical  difficulties  involved  with 

estimating  parameters  for  PH  distribution  and  the  fact  that  the  family  of  PH  distributions 

i. 
is  so  large  it  is  not  obvious  which  procedure  is  best.  The  criterion  for  comparison  is  the 

actual  average  cost  per  unit  time,  AACN,  after  N  replacements. 
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B.      COMPARISON  OF  NONPARAMETRIC  AND  PARAMETRIC  PROCEDURE 
FOR  THE  PH  DISTRIBUTION 

Simulation  was  used  to  compare  sequential  estimation  based  on  PH  distributions 
with  nonparametric  sequential  estimation,  Q  =  100,  C2  =  500.  System  lifetimes  were 
generated  by  an  increasing  failure  rate  PH  distribution  with  representation  (T41,a41) 
where, 


r41  = 


-5     4      1 
1-6     5 

11-7 


aA1  =    (1.0,0.0,0.0) 


For  this  representation  the  long  run  expected  cost  per  unit  time  is  given  in  Figure  3. 

To  give  the  parametric  procedure  a  chance,  the  first  25  system  lifetimes  are 
uncensored.  After  the  first  25  observations,  the  parameters  of  the  PH  distribution  are 
estimated  sequentially  as  in  Chapter  in  Section  B.  The  actual  average  costs  (AAC)  are 
calculated  and  compared  for  small,  medium  and  large  sample  sizes  with  planned  cost  Q 
=  100  and  unplanned  cost  C2  =  500.  The  20  initial  data  sets  are  simulated  as  shown  in 
the  Appendix  A. 

The  programs  are  separated  into  5  parts.  The  first  one  is  written  in  GAMS  as 
shown  in  Appendix  A,  to  estimate  the  parameters  of  a  PH  distribution  by  the  nonlinear 
programming  approach.  The  rest  are  written  in  Fortran  as  shown  in  Appendix  B,  to 
estimate  the  optimal  age  of  replacement  based  on  the  estimated  PH  distribution  from  the 
GAMS  program;  to  simulate  the  system  lifetime;  to  prepare  the  right  censored  data;  to 
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estimate  the  expected  lifetime,  second  and  third  standardized  moments;  and  finally  to 
tabulate  the  results  and  to  prepare  the  initial  data  for  the  next  estimation. 


LONG  RUN  EXPECTED  COST  PER  UNIT  TIME  OF  PH  DISTN. 

WITH  REPRESENTATION  T41,  a41  AND  C1  =  100,  C2=500 

o 

o 

CO 
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750 

EXPECTED 

700 
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z 
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LONG 
650 
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I                 i                 i                 i                 i                 i                 i                 i                 i 

10 

0                                       12                                      3                                      4 

REPLACEMENT  AGE 

Figure  3  Long  run  expected  cost  per  unit  time  curve  of  PH  distribution  with 
representation  T4,  and  <*41 

In  this  simulation,  the  number  of  transient  states  is  fixed  at  3  and  the  initial 

probability  vector  a  =  (1,  0,  0).  By  fixing  a  =  (1,  0,  0)  we  reduce  the  number  of 

parameters  to  be  estimated  and  use  a  model  in  which  all  systems  start  in  the  same  state, 

when  states  are  thought  of  as  the  level  of  system  repair,  this  choice  better  reflects  the 

idea  that  replacement  systems  are  as  good  as  new.  Bounded  parameters  and  user 
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interaction  are  used  in  guiding  each  estimation  during  the  simulation  to  a  preferable 
solution. 

Tables  1,  2  and  3  summarize  the  simulation  results  of  the  actual  average  cost 
resulting  from  sequential  estimation  of  PH  and  the  nonparametric  procedure.  Also  given 
are  the  signed  ranks  of  the  differences  in  the  actual  average  costs.  These  are  used  in  the 
Wilcoxon  Signed -Rank  test  to  test  the  hypothesis: 
H0:  E[AACP„]  =  E[AACNONP] 
Ha:  E[AACPH]  <  E[AACN0NP] 
by  using  T+  the  sum  of  the  positive  ranks  as  the  test  statistic  . 

Tables  1,  2  and  3  give  the  statistic  T+  as  35,  10  and  2  respectively.  When  compared  to 
the  value  in  Table  9  [Ref.  13:  pp.  780]  at  N=20  all  p-values  are  less  than  0.005  which 
implies  that  the  sequential  estimation  of  a  PH  distribution  is  better  than  sequential 
estimation  of  nonparametric  for  all  sample  sizes  (small,  medium  and  large).  The  box 
plots  in  Figures  5  and  6  of  actual  average  cost  versus  the  number  of  replacements  at 
different  sample  sizes  show  that  the  actual  average  costs  decrease  as  the  number  of 
replacements  increase  and  that  AACPH  for  the  parametric  procedure  decrease  more  than 
the  AACNONp  for  the  nonparametric  procedure. 

A  second,  PH  distribution  was  chosen  to  compare  the  parametric  PH  procedure 
with  the  nonparametric  procedure.  The  second  PH  distribution  has  an  average  long  run 
expected  cost  function  that  is  more  shallow  than  the  first  PH  distribution.  It  has 
representation  T42  and  a42 
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where 


T     = 


-5  4  1 
1-6  5 
12-7 


a..  =    (1.0,0.0,0.0)      , 


and  the  long  run  expected  cost  per  unit  time  is  given  in  Figure  4. 


LONG  RUN  EXPECTED  COST  PER  UNIT  TIME  OF  PH  DISTN. 
WITH  REPRESENTATION  T42.  or42  AND  01=100.  02  =  500 
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Figure  4  Long  run  expected  cost  per  unit  time  curve  of  PH  distribution  with 
representation  T^  and  att 

Tables  4,  5  and  6  summarize  the  result  and  give  the  statistic  T+  of  93,  81  and  47 
respectively.  The  p- value  taken  from  table  9  [Ref.  13:  pp.  780]  indicated  that  the 
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parametric  procedure  is  not  better  than  the  nonparametric  procedure  for  small  and 
moderate  sample  sizes.  In  fact,  for  some  cases  (e.g.  runs  16  and  20  in  Table  4  and  runs 
1,  5  and  16  in  Table  5)  the  actual  average  cost  resulting  from  using  the  nonparametric 
procedure  can  be  quite  a  bit  less  than  the  actual  average  cost  resulting  from  the 
parametric  procedure.  This  is  due  to  the  fact  that  C(</>)  is  shallow  at  <f>*,  thus  even  though 
C(4>)  may  be  a  reasonable  estimator  of  C(<f>),  the  variance  of  its  minimizing  value  $*  is 
higher  than  for  the  previous  PH  distribution.  This  means  that  </>*  is  often  underestimated. 
Because  C(</>)  increases  so  rapidly  to  the  left  of  <f>*,  underestimating  <j>*  can  increase  the 
actual  average  costs  dramatically.  Another  difficulty  is  that  the  estimated  C(<f>)  may  be 
relatively  flat  around  </>*.  This  causes  numerical  difficulties  with  the  nonlinear 
programming  algorithm.  For  large  sample  sizes  the  parametric  procedure  does  do  better 

A 

than  the  nonparametric  procedure.  Here  C(<f>)  estimates  C(<f>)  more  closely  and  it  is  easier 
for  the  nonlinear  programming  algorithm  to  identify  the  correct  solution. 

In  Figures  5  and  6,  the  improvement  in  actual  average  cost  with  sample  size  is 
shown  for  both  the  parametric  and  nonparametric  procedure  for  the  PH  distribution  with 
representation  T41  and  a41.  As  expected,  both  procedures  improve  (actual  average  costs 
decrease)  with  sample  size.  However,  the  actual  average  cost  for  the  parametric 
procedure  decreases  faster  than  for  the  nonparametric  procedure.  Both  procedures  exhibit 
considerable  variability.  The  solutions  have  got  still  mixed  that  we  should  do  more  study 
on  another  representation  until  the  more  precise  solutions  come  up. 
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Table  1  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  SMALL  SAMPLE 
SIZES  (N  =  50)  OF  PH  DISTRIBUTION  WITH  T41  AND  a41 


Run 

./\/\  v^pjj 

AACNONP 

/VrW^pp 

DIFFERENCE 

RANK 

1 

690.095 

754.555 

759.142 

-64.46 

15 

2 

714.728 

775.712 

780.872 

-60.444 

13 

3 

700.342 

684.379 

688.578 

15.963 

4  * 

4 

666.014 

725.331 

738.987 

-59.317 

12 

5 

714.330 

757.528 

735.058  $ 

-43.198 

8 

6 

544.322 

621.505 

621.507 

-77.183 

19 

7 

721.505 

810.149 

810.149 

-88.644 

20 

8 

565.854 

617.458 

636.556 

-51.604 

11 

9 

676.511 

725.394 

702.676  $ 

-48.883 

10 

10 

698.572 

766.182 

809.629 

-67.61 

17 

11 

814.989 

767.586 

801.150 

47.403 

9 

12 

644.921 

667.553 

672.116 

-22.632 

5 

13 

761.516 

686.794 

679.892  ! 

74.722 

18  * 

14 

608.477 

601.215 

617.711 

7.262 

1    * 

15 

706.591 

767.782 

782.198 

-61.191 

14 

16 

629.294 

667.805 

682.719 

-38.511 

6 

17 

688.572 

755.475 

748.353  $ 

-66.903 

16 

18 

636.963 

624.132 

665.344 

12.831 

3   * 

19 

583.226 

625.161 

609.036  $ 

-41.935 

7 

20 

644.065 

654.691 

676.005 

-10.626 

2 

J 

VACop  =  Actual 

average  cost  n 

jplacing  when  i 

ailure  occurs 

DIFFERENCE  =  AACPH  -  AAC 


'NONP 


*  positive  rank 

T+  =  35  ; 

p-value  <  0.005 
$  AACrj,  is  less  than  only  AACNONP 
!  AACrp  is  less  than  AACPH  and  AACNONP 
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Table  2   COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  MEDIUM 

SAMPLE  SIZES  (N  =  100)  OF  PH  DISTRIBUTION  WITH  T41  AND  a41 


Run 

J\J\\~*pii 

AACNONP 

.rv/W^pp 

DIFFERENCE 

RANK 

1 

648.127 

723.170 

732.263 

-75.043 

18 

2 

690.971 

786.554 

789.767 

-95.583 

20 

3 

649.886 

644.710 

683.933 

5.176 

2  * 

4 

696.858 

747.640 

755.982 

-50.782 

11 

5 

690.195 

688.749 

682.840  ! 

1.446 

1  * 

6 

577.742 

607.414 

610.850 

-29.672 

6 

7 

686.414 

766.365 

762.245  $ 

-79.951 

19 

8 

550.106 

612.046 

638.286 

-61.94 

14 

9 

637.268 

696.249 

689.552  $ 

-58.981 

12 

10 

684.253 

730.138 

786.920 

-45.885 

10 

11 

802.039 

790.698 

807.530 

11.341 

4  * 

12 

615.087 

687.719 

695.174 

-72.632 

17 

13 

677.302 

666.119 

674.830 

11.183 

3 

14 

586.838 

628.687 

638.451 

-41.849 

8 

15 

726.213 

749.507 

737.792  $ 

-23.294 

5 

16 

626.605 

694.232 

667.462  $ 

-67.627 

16 

17 

671.908 

732.053 

758.912 

-60.145 

13 

18 

632.954 

667.156 

727.346 

-34.202 

7 

19 

635.989 

679.399 

667.141  $ 

-43.41 

9 

20 

590.881 

657.064 

701.735 

-66.183 

15 

J 

VACRp  =  Actual 

average  cost  rei 

placing  when 

failure  occurs 

DIFFERENCE  -  AACPH  -  AACNONP 
*  positive  rank 

T+  =  10   ; 

p-value  <  0.005 
$  AACrf  is  less  than  only  AACNOnp 
!  AACrf  is  less  than  AACPH  and  AACNONP 
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Table  3  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  LARGE  SAMPLE 
SIZES  (N  =  200)  OF  PH  DISTRIBUTION  WITH  T41  AND  a4l 


Run 

l\J\\^-pyi 

AACNONP 

/Vr\V-^jyj 

DIFFERENCE 

RANK 

1 

622.560 

655.891 

655.865  $ 

-33.331 

7 

2 

658.640 

711.831 

708.966  $ 

-53.191 

11 

3 

623.823 

635.065 

671.611 

-11.242 

1 

4 

678.141 

749.866 

754.490 

-71.725 

16 

5 

674.548 

714.291 

707.926  $ 

-39.743 

9 

6 

585.435 

645.920 

647.927 

-60.485 

13 

7 

647.302 

729.159 

733.860 

-81.857 

20 

8 

581.711 

607.068 

649.386 

-25.357 

5 

9 

638.465 

678.078 

719.904 

-81.439 

19 

10 

655.709 

717.569 

766.020 

-61.86 

14 

11 

760.478 

779.362 

762.986  $ 

-18.884 

3 

12 

589.088 

649.545 

681.004 

-60.457 

12 

13 

696.081 

679.816 

712.761 

16.265 

2  * 

14 

602.032 

643.703 

649.832 

-47.8 

10 

15 

685.977 

767.133 

755.331  $ 

-81.156 

18 

16 

646.752 

723.719 

663.489  $ 

-76.967 

17 

17 

680.457 

713.146 

754.616 

-32.689 

6 

18 

674.300 

697.304 

740.394 

-23.004 

4 

19 

636.649 

674.845 

670.146$ 

-38.196 

8 

20 

588.468 

654.994 

670.194 

-66.526 

15 

J 

VACrf  =  Actui 

d  average  cost  i 

eplacing  when  f 

ailure  occurs 

DIFFERENCE  =  AACPH  -  AACNONP 
*  positive  rank 

T+  =  2   ; 

p-value  <  0.005 
$  AACrp  is  less  than  only  AACNOnp 
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Table  4  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  SMALL  SAMPLE 
SIZES  (N  =  50)  OF  PH  DISTRIBUTION  WITH  T^  AND  a^ 


Run 

A^\v^ppj 

AACNONP 

./V/W^-jvc 

DIFFERENCE 

RANK 

1 

618.463 

561.050 

545.871  ! 

57.413 

16* 

2 

604.098 

684.564 

686.650 

-80.466 

19 

3 

674.353 

664.878 

637.614  ! 

9.475 

6* 

4 

605.528 

559.966 

575.315 

45.562 

15  * 

5 

580.054 

521.136 

535.683 

58.918 

17  * 

6 

547.534 

553.458 

556.857 

-5.924 

4 

7 

632.508 

625.711 

625.711 

6.797 

5  * 

8 

632.510 

647.089 

651.264 

-14.579 

7 

9 

490.099 

578.531 

566.695  $ 

-88.432 

20 

10 

629.150 

668.268 

670.786 

-39.118 

13 

11 

535.032 

536.771 

535.194 

-1.739 

3 

12 

605.358 

644.458 

645.689 

-39.100 

12 

13 

|    519.537 

552.440 

552.443 

-32.903 

10 

14 

532.469 

553.664 

569.868 

-21.195 

8 

15 

671.880 

672.507 

672.160$ 

-0.627 

1 

16 

645.125 

571.650 

576.024 

73.475 

18  * 

17 

495.932 

533.334 

539.149 

-37.402 

11 

18 

648.741 

677.909 

682.898 

-29.168 

9 

19 

577.580 

576.101 

600.152 

1.479 

2  * 

20 

560.629 

519.418 

532.841 

41.211 

14  * 

AACrf  =  Actual  average  cost  replacing  when  failure  occurs 

DIFFERENCE  =  AACPH  -  AACNONp 
*  positive  rank 

T+  =  93  ; 

p-value  >  0.05 
$  AACrf  is  less  than  only  AACNONP 
!  AAC^  is  less  than  AACPH  and  AACNONP 
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Table  5   COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  MEDIUM 

SAMPLE  SIZES  (N  =  100)  OF  PH  DISTRIBUTION  WITH  T^  AND  a^ 


Run 

f\J\\^T>lI 

AACN0NP 

/vr\^pp 

DIFFERENCE 

RANK 

1 

573.310 

548.751 

539.868  ! 

24.559 

13  * 

2 

588.501 

633.798 

623.661  $ 

-45.297 

16 

3 

618.891 

626.148 

608.227  ! 

-7.257 

2 

4 

621.602 

590.532 

601.378 

31.070 

14  * 

5 

557.189 

521.040 

535.590 

36.149 

15  * 

6 

549.717 

562.877 

569.251 

-13.160 

5 

7 

651.656 

629.631 

630.823 

22.025 

10* 

8 

615.798 

624.193 

628.347 

-8.395 

4 

9 

475.757 

574.707 

567.688  $ 

-98.95 

20 

10 

602.762 

626.458 

628.035 

-23.696 

12 

11 

535.426 

537.073 

549.888 

-1.647 

1 

12 

580.480 

598.790 

590.995  $ 

-18.310 

7 

13 

525.770 

597.208 

597.210 

-71.438 

19 

14 

518.669 

540.978 

553.372 

-22.309 

11 

15 

634.978 

627.571 

633.019 

7.407 

3  * 

16 

675.027 

616.239 

618.041 

58.788 

18  * 

17 

535.911 

585.289 

587.998 

-49.378 

17 

18 

617.194 

636.912 

645.975 

-19.718 

9 

19 

550.468 

564.214 

579.831 

-13.746 

6 

20 

539.256 

520.034 

531.189 

19.222 

8  * 

J 

VACrf  =  Acta 

d  average  cost  i 

eplacing  when 

failure  occurs 

DIFFERENCE  =  AACPH  -  AACNONP 
*  positive  rank 

T+  =  81    ; 

p-value  >  0.05 
$  AACrf  is  less  than  only  AACNONP 
!  AACrj,  is  less  than  AACPH  and  AACNQNP 
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Table  6  COMPARISON  OF  ACTUAL  AVERAGE  COSTS  FOR  LARGE  SAMPLE 
SIZES  (N  =  200)  OF  PH  DISTRIBUTION  WITH  T^  AND  a^ 


Run 

/VrVv^pjj 

AACNONP 

/VfW^pp 

DIFFERENCE 

RANK 

1 

556.373 

550.131 

544.302  ! 

6.242 

3  * 

2 

574.206 

616.553 

621.391 

-42.347 

17 

;   3 

578.172 

591.186 

606.516 

-13.014 

7 

4 

594.573 

579.915 

588.547 

14.658 

8  * 

!     5 

!    580.261 

551.976 

560.463 

28.285 

14* 

6 

526.351 

558.604 

573.294 

-32.253 

15 

7      1 

632.238 

622.079 

624.275 

10.159 

6* 

8 

612.338 

621.511 

610.932  ! 

-9.173 

5 

9 

460.437 

567.069 

563.305  $ 

-106.632 

20 

10 

557.424 

599.181 

602.122 

-41.757 

16 

11 

555.711 

539.303 

551.773 

16.408 

9  * 

12 

560.534 

584.128 

589.078 

-23.594 

12 

13 

576.939 

621.856 

621.857 

-44.917 

18 

14 

568.925 

595.619 

599.305 

-26.694 

13 

15 

562.097 

583.057 

595.668 

-20.960 

11 

16 

584.778 

604.211 

606.979 

-19.433 

10 

17 

549.979 

607.061 

610.878 

-57.082 

19 

18 

631.719 

628.295 

639.159 

3.424 

2  * 

19 

578.172 

577.936 

607.344 

0.236 

1  * 

20 

538.683 

530.002 

570.516 

8.681 

4  * 

i 

^ACnp  =  Actuj 

il  average  cost  i 

eplacing  when 

failure  occurs 

DIFFERENCE  =  AACPH  -  AACNONP 
*  positive  rank 

T+  =  47   ; 

0.01  <  p-value  <  0.025 
$  AACrf  is  less  than  only  AACNONP 
!  AACrf  is  less  than  AACPH  and  AACNONP 
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Figure  5  The  box  plot  of  actual  average  cost  for  different  sample  sizes  of 
phase  type  sequential  procedure  for  T4,  and  a41 
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Figure  6  The  box  plot  of  actual  average  cost  for  different  sample  sizes  of 
nonparametric  sequential  procedure  for  T41  and  a4, 
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V.CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  the  age  replacement  policy  has  been  considered.  A  system  is  replaced 
at  the  time  of  failure  or  at  a  scheduled  replacement  age  whichever  comes  first.  The 
replacement  system  is  assumed  to  be  as  good  as  new.  The  objective  is  to  achieve  the 
minimum  long  run  expected  maintenance  cost  per  unit  time.  The  system  lifetimes  are 
assumed  to  have  PH  distributions.  Estimating  the  parameters  of  a  PH  distribution  is 
difficult  because  it  may  have  many  parameters  and  its  density  function  involves  a  matrix 
exponential.  Moment  matching  with  a  nonlinear  programming  approach  is  used  for 
estimating  the  parameters.  This  method  does  not  guarantee  that  the  estimated  cost 
function  will  have  a  unique  and  finite  minimum.  In  this  thesis  we  restric  our  attention  to 
the  case  with  initial  probability  vector  a  =  (1,  0,  0).  Sequential  estimation  using  the 
phase  type  procedure  gives  smaller  average  costs  than  the  nonparametric  procedure  for 
both  small  and  large  sample  sizes  when  the  underlying  long  ran  average  cost  function  has 
a  very  distinct  minimum  at  <f>*.  When  this  cost  function  is  shallow  around  <f>*,  the 
parametric  procedure  does  not  do  better  than  the  nonparametric  procedure  for  small 
samples.  In  fact,  for  cost  functions  that  are  shallow,  it  is  better  not  to  use  a  maintenance 
policy  i.e.  to  take  $*  =  oo  t  until  the  sample  sizes  are  large  enough  to  give  reasonable 
estimates  of  the  underlying  parameters.  The  reason  for  this  is  that  for  such  cost 
functions,  overestimating  <f>*  does  not  increase  the  long  ran  average  cost  much,  at  the 
same  time  the  decrease  in  the  amount  of  censoring  with  over  estimating  <i>*  greatly 
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improves  estimation  of  the  underlying  distribution  function  F.  Conversely 
underestimating  <f>*  for  these  cost  functions  causes  a  drastic  increase  in  long  run  average 
costs  and  increase  censoring  making  estimation  very  difficult. 

Recommendation  for  the  future  research 
The  following  forms  a  list  of  future  research  initiatives: 

-Examine  the  impact  of  the  product-limit  estimator  in  estimating  standardized 
moments  on  the  rest  of  procedure. 

-Do  more  experimentation  on  different  representations  of  phase  type  distributions 
to  get  more  precise  conclusion. 

-Look  for  other  policies  using  phase-sensitive  estimate  and  modified  replacement 
costs  over  time. 

-The  phase-type  sequential  estimations  we  have  been  so  far  assume  the  underlying 
lifetimes  are  iid  and  after  replacement  the  system  is  new.  To  be  more  realistic,  we  can 
modify  the  initial  probability  vector,  increase  transition  rates  between  states,  or  both, 
over  time. 

-Use  phase-type  sequential  estimation  when  the  underlying  life  distribution  F  comes 
from  Gamma,  Weibull  etc.,  that  have  increasing  failure  rate  and  compare  with  the 
nonparametric  procedure.  Since  we  can  use  the  denseness  of  phase  type  distribution 
property  to  approximate  the  set  of  distribution  with  support  on  [0,  a)  i.e.,  the  set  of 
lifetimes,  we  may  get  a  better  method  when  the  distribution  of  lifetimes  is  not  known. 
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APPENDIX  A. 

$TTTLE       PH  DISTRIBUTION  ESTIMATION  PARAMETERS 

* 

* GAMS  AND  DOLLAR  CONTROL  OPTIONS 

*  (See  Appendices  B  &  C) 

OPTIONS 

WORK  =  100000, 

LIMCOL  =  0  ,  LJMROW    =  0  ,  SOLPRINT  =  OFF  ,  DECIMALS  =  2 

RESLIM  =  100,  ITERLIM  =100000,  OPTCR  =  0.0  ,  SEED    =  3141; 


* DEFINITIONS  AND  DATA 

*  This  program  uses  nonlinear  programming  to  estimate  the  3-transient  states  of  phase 

*  type  distribution  parameters.  Matching  the  second  and  third  standardized  moments  are 

*  used.  All  adjustable  data  are  prepared  by  the  "FILE  SCALAR2"  and  uses  the 

*  "$INCLUDE"  statement  bring  to  the  program. 

$INCLUDE  "FILE  SCALAR2" 


* -MODEL- 

VARIABLE 

Z  objective  function  value 

A  first  entry  of  the  matrix 

B  second  entry  of  the  matrix 

C  third  entry  of  the  matrix 

D  fourth  entry  of  the  matrix 

E  fifth  entry  of  the  matrix 

F  sixth  entry  of  the  matrix 

G  seventh  entry  of  the  matrix 

H  eighth  entry  of  the  matrix 

I  ninth  entry  of  the  matrix 

AL1  initial  probability  of  being  in  state  1 

AL2  initial  probability  of  being  in  state  2 

AL3  initial  probability  of  being  in  state  3  ; 

EQUATION 

OBJ  Define  objective  function 

EQ1   Set  initial  probability  of  being  in  state  l(nonnegative) 

EQ2  Set  initial  probability  of  being  in  state  2(nonnegative) 
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EQ3  Set  initial  probability  of  being  in  state  3(nonnegative) 

EQ4  Sum  over  all  of  initial  prob.  less  than  or  equal  1.0 

EQ5   Set  the  first  diagonal  entry  to  be  nonpositive 

EQ6  Set  the  second  diagonal  entry  to  be  nonpositive 

EQ7  Set  the  third  diagonal  entry  to  be  nonpositive 

EQ8   Set  off  diagonal  entries  in  the  first  row  to  be  positive 

EQ9  Set  off  diagonal  entries  in  the  first  row  to  be  positive 

EQ10  Set  sum  off  diagonal  of  the  first  row  less  than  or  equal  first  diagonal 

EQ11  Set  off  diagonal  entries  in  the  second  row  to  be  positive 

EQ12  Set  off  diagonal  entries  in  the  second  row  to  be  positive 

EQ13  Set  sum  off  diagonal  of  the  second  row  less  than  or  equal  second  diagonal 

EQ14  Set  off  diagonal  entries  in  the  third  row  to  be  positive 

EQ15  Set  off  diagonal  entries  in  the  third  row  to  be  positive 

EQ16  Set  sum  off  diagonal  of  the  third  row  less  than  or  equal  third  diagonal 

EQ17  Define  the  expected  lifetime   ; 

♦MINIMIZE 
OBJ.. 
Z  =E= 

***   3*SQR(C  -  C(0))   *** 

3*POWER((((2*ALl*(POWER((-(F*H)+E*I),2)+(C*H-B*I)*(F*G-D*I) 

+ (-(C*E)  +B*F)*(-(E*G) +D*H) 
+ (-(F*H) +E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) + (-(C*E) +B*F)*(B*G-A*H) 
+  (-(F*H)+E*I)*(-(C*E)+B*F)+(C*H-B*I)*(C*D-A*F)  +  (-(C*E)+B*F)*(-(B*D)  + 
A*E)) 

+  2*AL2*((F*G-D*I)*(-(F*H)+E*I)+(-(C*G)+A*I)*(F*G-D*I)+(-(E*G)+D*H)* 

(C*D-A*F) 
+  (F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H) 
+ (F*G-D*I)*(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)  *(-(B*D) + A*E)) 

+  2*AL3*((-(E*G)+D*H)*(-(F*H)+E*I)+(B*G-A*H)*(F*G-D*I)  +  (-(B*D)+A*E)* 

(-(E*G)+D*H) 
+  (-(E*G) +D*H)*(C*H-B*I)  +  (B*G-A*H)*(-(C*G) + A*I) + (-(B*D) + A*E)* 
(B*G-A*H) + (-(E*G) +D*H)*(-(C*E)  +B*F)  +  (B*G-A*H)*(C*D-A*F)  + 
POWER((-(B*D) + A*E),2)))/ 

POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),2) 

-POWER((-(ALl*(-F*H+E*H-C*H-B*I-(C*E)+B*F) 

+AL2*(F*G-D*I-C*G+  A*I+C*D-A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 
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(-(C*E*G)+B*F*G+C*D*H-A*F*H-B*D^+A*E*I)),2))**0.5/ 

(-(AL1*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F) 

+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 

+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 

(-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I))  -  MM2),2) 

***  W2*SQR(GAMMA  -  GAMMA(O))   *** 

+  POWER(((((-6*ALl*((POWER((-(F*H)+E*I),2)+(C*H-B*I)*(F*G-D*I)  + 
(-(C*E)+B*F)*(-(E*G)+D*H))*(-(F*H)+E*I) 

+((-(F*H) +E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I) +(-(C*E)  +B*F)*   (B*G-A*H))* 
(F*G-D*I) 

+ ((-(F*H) +E*I)*(-(C*E) +B*F) + (C*H-B*I)*(C*D-A*F)  +  (-(C*E) +B*F)* 
(-(B*D) + A*E))*(-(E*G) +D*H) 

+(POWER((-(F*H)+E*I),2)+(C*H-B*I)*(F*G-D*I)+(-(C*E)+B*F)'K  (-(E*G)+D*H))* 
(C*H-B*I) 

+ ((-(F*H) + E*I)  *(C*H-B*I) + (C*H-B*I)  ^-(C*G) + A*I) + (-(C*E)  +B*F)*(B*G-A*H))* 
(-(C*G)+A*I) 

+ ((-(F*H) +E*I)*(-(C*E) +B*F)  +  (C*H-B*I)*(C*D-A*F) + (-(C*E) +B*F)* 
(-(B*D) + A*E))*(B*G-A*H) 

+  (POWER((-(F*H) +E*I),2)  +  (C*H-B*I)*(F*G-D*I) + (-(C*E) +B*F)* 
(-(E*G) +D*H))*(-(C*E) +B*F) 

+((-(F*H) +E*I)*(C*H-B*I)  +  (C*H-B*I)*(-(C*G) + A*I) +(-(C*E)  +B*F)* 
(B*G-A*H))*(C*D-A*F) 

+((-(F*H) +E*I)*(-(C*E) +B*F) + (C*H-B*I)*(C*D-A*F) + (-(C*E) +B*F)* 
(-(B*D) + A*E))*(-(B*D) + A*E)) 

-6*AL2*(((F*G-D*I)*(-(F*H) +E*I) + (-(C*G)  4-  A*I)  *(F*G-D*I) + (-(E*G)  +D*H)* 
(C*D-A*F))*(-(F*H) +E*I) 

+((F*G-D*I)*(C*H-B*I)+POWER((-(C1,cG)+A*I),2)+(C*D-A*F)*(B*G-A*H))* 
(F*G-D*I) 

+ ((F*G-D*I)*(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)* 
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(-(B*D) + A*E))*(-(E*G) +D*H) 

+ ((F*G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I)*(F*G-D*I) + (-(E*G) +D*H)* 
(C*D-A*F))*(C*H-B*I) 

+((F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H))* 
(-(C*G)+A*I) 

+((F*G-D*I),,e(-(C*E)+B*F)+(-(C*G)+A*I)*(C*D-A*F)+(C*D-A*F)* 
(-(B*D) + A*E))*(B*G-A*H) 

+  ((F*G-D*I)*(-(F*H) +E*I) + (-(C*G) + A*I)*(F*G-D*I)  +  (-(E*G) +D*H)* 
(C*D-A*F))*(-(C*E)+B*F) 

+((F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H))* 
(C*D-A*F) 

+ ((F*G-D*I)*(-(C*E)  +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)* 
(-(B*D) + A*E))*(-(B*D) + A*E)) 

-6*AL3*(((-(E*G) +D*H)*(-(F*H) +E*I)  4-  (B*G-A*H)*(F*G-D*I) + (-(B*D) + A*E)* 
(-(E*G) +D*H))*(-(F*H) +E*I) 

+ ((-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*I)  +  (-(B*D) + A*E)* 
(B*G-A*H))*(F*G-D*I) 

+  ((-(E*G) +D*H)*(-(C*E) +B*F) + (B*G-A*H)*(C*D-A*F)  + 
POWER((-(B*D) + A*E),2))*(-(E*G) +D*H) 

+ ((-(E*G) +D*H)*(-(F*H) +E*I)  +  (B*G-A*H)*(F*G-D*I) + (-(B*D) + A*E)* 
(-(E*G)+D*H))*(C*H-B*I) 

+ ((-(E*G) +D*H)*(C*H-B*I) + (B*G-A*H)*(-(C*G) + A*I)  +  (-(B*D) + A*E)* 
(B*G-A*H))*(-(C*G)+A*I) 

+  ((-(E*G) +D*H)*(-(C*E)  +B*F) + (B*G-A*H)*(C*D-A*F)  + 
POWER((-(B*D)+A*E),2))*(B*G-A*H) 

+ ((-(E*G) +D*H)*(-(F*H) +E*I)  +  (B*G-A*H)*(F*G-D*I)  4-  (-(B*D) + A*E)* 
(-(E*G)+D*H))*(-(C*E)+B*F) 

+((-(E*G)+D*H)*(C*H-B*I)+(B*G-A*H)*(-(C*G)+A*I)  +  (-(B*D)+A*E)* 
(B*G-A*H))*(C*D-A*F) 
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+ ((-(E*G) +D*H)*(-(C*E) +B*F) + (B*G-A*H)*(C*D-A*F)  + 
POWER((-(B*D) + A*E)  ,2))  *(-(B*D) + A*E))) 

/POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),3) 

+3*((AL1*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F)+AL2*(F*G-D*I-C*G+A*I+C*D- 
A*F)+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/(-C*E*G+B*F*G+C*D*H-A*F 
*H-B*D*I+A*E*I))* 

***  expected  of  lifetime  square   *** 

((2*AL1  *(POWER((-(F*H) +E*I),2)  +  (C*H-B*I)*(F*G-D*I) + (-(C*E) +B*F)* 
(-(E*G) +D*H) + (-(F*H) +E*I)*(C*H-B*I) + (C*H-B*I)*(-(C*G) + A*I)  + 
(-(C*E)+B*F)*(B*G-A*H)+(-(F*H)+E*I)*(-(C*E)+B*F)  +  (C*H-B*I)*(C',eD-A*F) 
+  (-(C*E) +B*F)*(-(B*D) + A*E)) 

+2*AL2*((F*G-D*I)*(-(F*H) +E*I)  +  (-(C*G) + A*I)*(F*G-D*I) 

+  (-(E*G)  +D*H)*(C*D-A*F)+(F*G-D*I)*(C*H-B*I)+POWER((-(C*G)+A*I),2) 

+ (C*D-A*F)*(B*G-A*H)  +  (F*G-D*I)*(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) 

+(C*D-A*F)*(-(B*D)+A*E)) 

+2*AL3*((-(E*G) +D*H)*(-(F*H) +E*I)  +  (B*G-A*H)*(F*G-D*I) 
+ (-(B*D) + A*E)*(-(E*G) +D*H) + (-(E*G) +D*H)*(C*H-B*I)  +  (B*G-A*H)* 
(-(C*G) + A*I) + (-(B*D) + A*E)*(B*G-A*H)  +  (-(E*G) +D*H)*(-(C*E) +B*F) 
+  (B*G-A*H)*(C*D-A*F) +POWER((-(B*D) + A*E),2)))/ 

POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),2)) 

***  2  time  cube  of  the  expected  lifetime  *** 

-2*POWER(((ALl*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F) 

-f-AL2*(F*G-D*I-C*G+A*I+C*D-A*F)+AL3*(-E,,eG+D*H+B*G-A*H-B*D+A*E))/ 
(-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I)),3)) 

/((2*  AL1  *(POWER((-(F*H) + E*I)  ,2) + (C  *H-B*I)  *(F*G-D*I) + (-(C*E)  +B*F)  * 
(-(E*G) +D*H)  +  (-(F*H) +E*I)*(C*H-B*I) +(C*H-B*I)*(-(C*G) + A*I)  + 
(-(C*E)+B*F)*(B*G-A*H)+(-(F*H)+E*I)*(-(C*E)+B*F)  +  (C*H-B*I)*(C*D-A*F) 
+(-(C*E)+B*F)*(-(B*D)+A*E)) 

+2*AL2*((F*G-D*I)*(-(F*H) +E*I) +(-(C*G) + A*I)*(F*G-D*I)  + 
(-(E*G) +D*H)*(C*D-A*F) + (F*G-D*I)*(C*H-B*I)  + 
POWER((-(C*G)+A*I),2)+(C*D-A*F)*(B*G-A*H)+(F*G-D*I)* 
(-(C*E) +B*F) + (-(C*G) + A*I)*(C*D-A*F) + (C*D-A*F)*(-(B*D) + A*E)) 
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+  2*AL3*((-(E*G) +D*H)*(-(F*H)  +E*I) + (B*G-A*H)*(F*G-D*I)  + 

(-(B*D)  +  A*E)*(-(E*G)  +  D*H)  +  (-(E*G)  +  D*H)*(C*H-B*I)  +  (B*G-A*H)* 

(-(C*G)  +  A*I)  +  (-(B*D)  +  A*E)*(B*G-A*H)  +  (-(E*G)  +  D*H)*(-(C*E)  +  B*F) 

+(B*G-A*H)*(C*D-A*F)+POWER((-(B*D)+A*E),2)))/ 

(POWER((-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I),2)) 

-P0WER(((AL1*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F) 
+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 
+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 
(-(C*E*G)  +  B*F*G  +  C*D*H  -  A*F*H  -  B*D*I  +  A*E*I)),2))**1.5)  -  MM3),2); 

*  SUBJECT  TO 

EQ1.. 

AL1  =G=  0  ; 
EQ2.. 

AL2  =G=  0  ; 


EQ3.. 

AL3  =G=  0  ; 

EQ4.. 

AL1  +  AL2  +  AL3  =L=  1  ; 

EQ5.. 
A 

=L=  0  ; 

EQ6.. 
E 

=L=  0   ; 

EQ7.. 
I 

=L=  0   ; 

EQ8.. 
B 

=G=  0  ; 

EQ9.. 
C 

=G=  0  ; 

EQ10.. 

B  +  C  =L=  -A  ; 

EQll.. 
D 

=G=  0; 

EQ12.. 
F 

=G=  0  ; 

EQ13.. 

D  +  F  =L=  -E  ; 

EQ14.. 
G 

=G=  0  ; 
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EQ15.. 

H     =G=  0 ; 
EQ16.. 

G  +  H  =L=  -I  ; 
EQ17.. 

-(AL1*(-(F*H)+E*I+C*H-B*I-(C*E)+B*F) 
+AL2*(F*G-D*I-C*G+A*I+C*D-A*F) 
+AL3*(-E*G+D*H+B*G-A*H-B*D+A*E))/ 
(-(C*E*G)+B*F*G+C*D*H-A*F*H-B*D*I+A*E*I)  =E=  EX  ; 


* REPORT- 

* 

MODEL  MATCHM  /ALL/  ; 

*  Parameters  bounding 

A.LO  =  -7.0;  B.LO  =  3.0;  C.LO  =  0.0;  D.LO  =  0.0;  E.LO  =-7.0; 
A.UP  =  -5.0;  B.UP  =  5.0;  CUP  =  2.0;  D.UP  =  1.0;  E.UP  =-5.0; 
F.LO  =   4.0;  G.LO  =  0.0;  H.LO  =  0.0;  I.LO  =  -8.0; 
F.UP  =   5.0;  G.UP  =  1.0;  H.UP  =  2.0;  I.UP  =  -7.0; 
AL1.LO  =  1.0;  AL2.LO  =  0.0;  AL3.LO  =  0.0; 
AL1.UP  =  1.0;  AL2.UP  =  0.0;  AL3.UP  =  0.0; 

*  The  initial  solution 

A.L  =  AA  ;  B.L  =  BB  ;  C.L  =  CC;  D.L  =  DD;  E.L  =  EE  ; 
F.L  =  FF  ;  G.L  =  GG  ;  H.L  =  HH  ;  I.L  =  II  ; 
AL1.L  =  ALP1  ;  AL2.L  =  ALP2  ;  AL3.L  =  ALP3  ; 

SOLVE  MATCHM  USING  NLP  MINIMIZING  Z  ; 
DISPLAY  Z.L,A.L  ,  B.L,  C.L,  D.L,  E.L,  F.L,  G.L,  H.L,  I.L, 
AL1.L,  AL2.L,  AL3.L; 

*  Write  output  to  the  CMS  file 

FILE  OUT1  /FILE  OBJ  A  /  ; 

PUT  OUT1  ; 

PUT  Z.L  ; 

FILE  OUT2  /FILE  ALPHA  A  /  ; 

PUTOUT2  ; 

PUT  A.L  ; 

FILE  OUT3  /FILE  BRAVO  A  /  ; 

PUTOUT3  ; 
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PUT  B.L  ; 

FILE  0UT4  /FILE  CHAREE  A  / 

PUT  OUT4 

PUT  C.L  ; 

FILE0UT5 

PUTOUT5 

PUT  D.L  ; 

FILEOUT6 

PUTOUT6 

PUT  E.L  ; 

FILEOUT7 

PUTOUT7 

PUT  F.L  ; 

FELEOUT8 

PUT0UT8 

PUT  G.L  ; 

FILE0UT9 

PUT0UT9 

PUT  H.L  ; 

FILE  OUT10  /FILE  INDIA  A  / 

PUT  OUT10  ; 

PUT  I.L  ; 

FILE  OUT11  /FILE  AL1  A  /  ; 

PUTOUT11  ; 

PUTAL1.L; 

FILE  OUT12  /FILE  AL2  A  /  ; 

PUT  OUT12  ; 

PUT  AL2.L  ; 

FILE  OUT13  /FILE  AL3  A  /  ; 

PUT  OUT13  ; 

PUT  AL3.L  : 


/FILE  DELTA  A  / 


/FILE  ECHO  A  / 


/FILE  FOXTROT  A  / 


/FILE  GOLF  A  / 


/FILE  HOTEL  A  / 
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THE  INITIAL  COMPLETE  SYSTEM  LIFETIME  SET  1-5 

N 

SET  1 

SET  2 

SET  3 

SET  4 

SET  5 

1 

0.1226 

0.0740 

0.0855 

0.0528 

0.1143 

2 

0.1455 

0.1413 

0.2157 

0.1139 

0.1999 

3 

0.2662 

0.2261 

0.2209 

0.2412 

0.2002 

4 

0.2925 

0.2629 

0.2839 

0.2892 

0.2422 

5 

0.3205 

0.2831 

0.2852 

0.3118 

0.3164 

6 

0.4214 

0.3340 

0.3116 

0.3485 

0.3398 

7 

0.4394 

0.3911 

0.3403 

0.3973 

0.3534 

8 

0.4685 

0.4617 

0.4517 

0.4442 

0.3893 

9 

0.5804 

0.4761 

0.4691 

0.4686 

0.4410 

10 

0.5871 

0.4909 

0.5254 

0.4777 

0.4766 

11 

0.6228 

0.5290 

0.5536 

0.4885 

0.4964 

12 

0.6534 

0.5818 

0.5993 

0.4970 

0.5074 

13 

0.7220 

0.5962 

0.6575 

0.5114 

0.5240 

14 

0.7590 

0.6233 

0.6623 

0.6073 

0.5825 

15 

0.7618 

0.6419 

0.6855 

0.6433 

0.6033 

16 

0.8206 

0.7407 

0.7156 

0.7576 

0.6255 

17 

0.8245 

0.7682 

0.7690 

0.8057 

0.7391 

18 

0.8317 

0.8312 

0.7866 

0.8740 

0.7723 

19 

1.0056 

0.8693 

0.8852 

1.0004 

0.9128 

20 

1.0259 

0.9727 

0.9473 

1.0424 

1.0832 

21 

1.0853 

1.1084 

1.0478 

1.0634 

1.1936 

22 

1.1659 

1.1272 

1.1155 

1.0659 

1.2003 

23 

1.1852 

1.2037 

1.2936 

1.2238 

1.8976 

24 

1.2040 

1.3319 

1.3906 

1.6404 

2.0769 

25 

1.3069 

2.5362 

1.6715 

3.0033 

2.2631 

E[X] 

0.705 

0.704 

0.679 

0.735 

0.742 

2nd  MM 

0.47917 

0.70936 

0.57558 

0.80566 

0.77290 

3rd  MM 

0.03190 

1.88276 

0.73612 

2.23985 

1.40912 
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THE  INE 

[TAL  COMPLI 

CTE  SYSTEM  LIFETIME  SET  6-10 

N 

SET   6 

SET    7 

SET    8 

SET   9 

SET   10 

1 

0.1337 

0.0402 

0.1759 

0.0830 

0.0761 

2 

0.1538 

0.0634 

0.2031 

0.1101 

0.0913 

3 

0.2452 

0.1939 

0.2106 

0.2083 

0.2235 

4 

0.2876 

0.1954 

0.2376 

0.2339 

0.3309 

5 

0.3591 

0.1958 

0.3043 

0.2408 

0.3470 

6 

0.4589 

0.2353 

0.3719 

0.2726 

0.3662 

7 

0.5161 

0.2775 

0.3850 

0.2777 

0.4091 

8 

0.6030 

0.3185 

0.3892 

0.2894 

0.4523 

9 

0.6181 

0.3579 

0.4032 

0.3137 

0.5581 

10 

0.6761 

0.3922 

0.4199 

0.3227 

0.5800 

11 

0.6851 

0.3985 

0.6029 

0.3390 

0.5805 

12 

0.7738 

0.5206 

0.6406 

0.3873 

0.5993 

13 

0.8541 

0.5233 

0.6506 

0.3878 

0.6354 

14 

0.8654 

0.6261 

0.7156 

0.4055 

0.6571 

15 

0.9539 

0.6336 

0.7584 

0.4076 

0.7529 

16 

1.0193 

0.6776 

0.7677 

0.4333 

0.7874 

17 

1.0352 

0.8705 

0.8754 

0.6178 

0.7991 

18 

1.0677 

1.1124 

0.8836 

0.6245 

0.8572 

19 

1.1196 

1.1228 

0.9575 

0.6455 

0.8920 

20 

1.1995 

1.1742 

0.9805 

0.9824 

0.9024 

21 

1.2301 

1.1872 

1.3367 

1.0054 

0.9699 

22 

1.3043 

1.1887 

1.4040 

1.1143 

0.9734 

23 

2.0273 

1.4849 

1.5844 

1.2279 

1.0452 

24 

3.2891 

1.7218 

1.6250 

1.2420 

1.1461 

25 

3.5360 

2.5301 

1.8392 

2.1319 

1.1594 

E[X] 

1.000 

0.722 

0.749 

0.572 

0.648 

2nd  MM 

0.82457 

0.81263 

0.63261 

0.80948 

0.46844 

3rd  MM 

1.89734 

1.25687 

0.79102 

1.69217 

0.14707 
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THE  INITIAL  COMPLETE  SYSTEM  LIFETIME  SET  11-15 

N 

SET    11 

SET   12 

SET   13 

SET   14 

SET   15 

1 

0.0268 

0.0433 

0.0470 

0.1730 

0.1101 

2 

0.0745 

0.1285 

0.0922 

0.2556 

0.1386 

3 

0.1235 

0.2052 

0.1295 

0.3124 

0.1709 

4 

0.1800 

0.2435 

0.2064 

0.3125 

0.1998 

5 

0.2526 

0.2924 

0.2924 

0.3447 

0.2075 

6 

0.2704 

0.3335 

0.2957 

0.3894 

0.2588 

7 

0.2972 

0.3645 

0.3110 

0.4044 

0.2669 

8 

0.3697 

0.4936 

0.3209 

0.5555 

0.2737 

9 

0.3811 

0.4984 

0.4173 

0.5796 

0.2905 

10 

0.3990 

0.5064 

0.4577 

0.5803 

0.3506 

11 

0.4070 

0.5778 

0.5168 

0.5928 

0.3604 

12 

0.4855 

0.6699 

0.5179 

0.6056 

0.4117 

13 

0.6402 

0.7442 

0.5258 

0.6525 

0.5809 

14 

0.7087 

0.7502 

0.5456 

0.7299 

0.5981 

15 

0.7408 

0.8656 

0.5547 

0.7355 

0.6558 

16 

0.8496 

0.8831 

0.5559 

0.8122 

0.7306 

17 

0.9290 

1.0607 

0.5633 

0.9379 

0.7388 

18 

0.9873 

1.1003 

0.5692 

0.9433 

0.7400 

19 

0.9983 

1.1834 

0.6501 

1.0197 

0.7786 

20 

1.0318 

1.2235 

0.7978 

1.2675 

0.7936 

21 

1.0937 

1.3370 

0.8750 

1.3842 

1.1532 

22 

1.1137 

1.5754 

1.0843 

1.6525 

1.2444 

23 

1.2706 

1.8085 

1.1683 

1.8276 

1.2781 

24 

1.4875 

1.9023 

1.1881 

1.8545 

1.3428 

25 

1.4958 

1.9894 

1.3242 

2.8764 

1.5937 

BEX] 

0.665 

0.831 

0.560 

0.872 

0.611 

2nd  MM 

0.64986 

0.66778 

0.60540 

0.71398 

0.68553 

3rd  MM 

0.32858 

0.60712 

0.69085 

1.51250 

0.77568 
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THE  INITIAL  COMPLETE  SYSTEM  LIFETIME  SET  16-20 

N 

SET   16 

SET   17 

SET   18 

SET   19 

SET  20 

1          ! 

0.1677 

0.0565 

0.0398 

0.1798 

0.0202 

2 

0.1819 

0.1811 

0.1402 

0.2037 

0.1920 

3 

0.2139 

0.1941 

0.1506 

0.2197 

0.2205 

4 

0.3436 

0.2810 

0.1792 

0.3159 

0.2209 

5 

0.3641 

0.3315 

0.1922 

0.3464 

0.3406 

6 

0.4393 

0.4147 

0.3787 

0.3629 

0.3645 

7 

0.4800 

0.4249 

0.4660 

0.3744 

0.3941 

8 

0.5328 

0.4629 

0.4944 

0.3793 

0.4130 

9 

0.5643 

0.4939 

0.5952 

0.3870 

0.5525 

10 

0.7220 

0.5382 

0.6198 

0.3875 

0.5729 

11 

0.7393 

0.5776 

0.6269 

0.4011 

0.5787 

12          | 

0.7666 

0.6003 

0.6276 

0.4906 

0.5835 

13           1 

0.7766 

0.6055 

0.6608 

0.7436 

0.6138 

14 

0.7964 

0.6230 

0.7183 

0.7690 

0.6389 

15 

0.7981 

0.6477 

0.7970 

0.7824 

0.6609 

16 

0.8228 

0.7631 

0.7997 

0.8226 

0.6624 

17 

0.8390 

0.7980 

0.8065 

0.8492 

0.7091 

18 

!      0.9422 

0.7988 

0.8065 

0.9542 

0.7178 

19 

1.0115 

0.8518 

0.9682 

1.2333 

0.7447 

20 

1.0381 

0.9091 

0.9939 

1.3601 

1.2027 

21 

1.1195 

1.0447 

1.0028 

1.4489 

1.3247 

22 

1.2716 

1.1973 

1.0875 

1.5989 

1.4412 

23 

1.3426 

1.2331 

1.2722 

2.2702 

1.6155 

24 

1.6082 

1.2339 

1.4913 

2.2971 

1.7162 

25 

1.7724 

1.3386 

1.8663 

3.1519 

2.1912 

EtX] 

0.786 

0.664 

0.711 

0.893 

0.748 

2nd  MM 

0.52328 

0.51958 

0.60066 

0.83695 

0.70314 

3rd  MM 

0.55452 

0.32242 

0.66254 

1.46868 

1.13073 
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APPENDIX  B. 

PROGRAM  COST 
C 

C  THIS  PROGRAM  PROVIDES  THE  OPTIMAL  AGE  REPLACEMENT  OF  THE 

C  PHASE  TYPE  DISTRIBUTION.  IT  WORKS  ONLY  SPECIFIC  CASE  OF  3 

C  PHASES  OF  TRANSIENT  STATES .  THE  PARAMETERS  ARE  AS  FOLLOWS : 
C         T(I,J)    =  ENTRIES  OF  THE  TRANSITION  MATRIX  OF  PH  DISTN. 
C         AL(I)    =  INITIAL  PROBABILITY  OF  BEING  IN  STATE  I 
C         CI  PLANNED  MAINTENANCE  COST 

C         C2  UNPLANNED  MAINTENANCE  COST 

C         L  MAXIMUM  LIFETIME  THAT  WANT  TO  CALCULATE 

C         DET     =  DETERMINANT  OF  TRANSITION  MATRIX 

C  SMALL  =  THE  OPTIMAL  AGE  REPLACEMENT 

C  CS       =  AGE  REPLACEMENT  COSTS 

* 

*   INITIAL  INPUT 

* 

INTEGER  I,J,K,N 

REAL  CI,  C2,  TI(3,3),  T(3,3),  DET,  AL(3),  L,  OBJ 

REAL  TS(3,3),  SMALL,  PI,  A,  IDEN(3,3),  S,  NUMER 

REAL  DENOM,  P(3,3),  Q(3,3),  CS(400),  EXPTS(3,3) 

DATA  C1,C2,L,TS/100.0,500.0,4.0,9*0.0  / 

N  =  0 

S  =  0.01 

OPEN(UNIT  =11, FILE  =  'ALPHA') 

OPEN(UNIT  =12,FILE  =  'BRAVO') 

OPEN(UNIT  =13,FILE  =  'CHAREE') 

OPEN(UNIT  =14,FILE  =  'DELTA') 

OPEN(UNIT  =15,FILE  =  'ECHO') 

OPEN(UNIT  =16,FILE  =  'FOXTROT') 

OPEN(UNIT  =17,FILE  =  'GOLF') 

OPEN(UNIT  =18,FDLE  =  'HOTEL') 

OPEN(UNTT  =19,FILE  =  'INDIA') 

OPEN(UNIT  =20,FDLE  =  'AL1') 

OPEN(UNIT  =21, FILE  =  'AL2') 

OPEN(UNIT  =22,FILE  =  'AL3') 

OPEN(UNIT  =23, FILE  =  'OBJ') 

READ(11,*)T(1,1) 

READ(12,*)T(1,2) 

READ(13,*)T(1,3) 
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READ(14,*)T(2,1) 

READ(15,*)T(2,2) 

READ(16,*)T(2,3) 

READ(17,*)T(3,1) 

READ(18,*)T(3,2) 

READ(19,*)T(3,3) 

READ(20,*)AL(1) 

READ(21,*)AL(2) 

READ(22,*)AL(3) 

READ(23,*)OBJ 

IF(OBJ  .GT.  0.05)  GOTO  555 

*  CALCULATION  COST  FOR  EACH  REPLACEMENT  TIME 

* 

DET  =  -(T(1,3)*T(2,2)*T(3,1)) 
/         +  T(1,2)*T(2,3)*T(3,1) 
/         +  T(1,3)*T(2,1)*T(3,2) 
/        -T(1,1)*T(2,3)*T(3,2) 
/        -T(1,2)*T(2,1)*T(3,3) 
/         +  T(1,1)*T(2,2)*T(3,3) 
PRINT*, 'DET  =  \DET 

11(1,1)  =  (-(T(2,3)*T(3,2))  +  T(2,2)*T(3,3))/DET 
11(1,2)  =  (T(1,3)*T(3,2)-T(1,2)*T(3,3))/DET 
TI(1,3)  =  (-(T(1,3)*T(2,2))  +  T(1,2)*T(2,3))/DET 
H(2,l)  =  (T(2,3)*T(3,1)-T(2,1)*T(3,3))/DET 
H(2,2)  =  (-(T(1,3)*T(3,1))  +  T(1,1)*T(3,3))/DET 
TI(2,3)  =  (T(1,3)*T(2,1)-T(1,1)*T(2,3))/DET 
11(3,1)  =  (-(T(2,2)*T(3,1))  +  T(2,1)*T(3,2))/DET 
H(3,2)  =  (T(1,2)*T(3,1)-T(1,1)*T(3,2))/DET 
TI(3,3)  =  (-(T(1,2)*T(2,1))  +  T(1,1)*T(2,2))/DET 
DO  101  =  1,3 
DO  5  J  =  1,3 

IF(I.EQ.J)THEN 

IDEN(I,J)  =  1.0 
ELSE 

IDEN(I,J)  =  0.0 
ENDEF 
5      CONTINUE 

10  CONTINUE 

* 

100  CONTINUE 
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N  =  N+l 
DO  20  I  =  1,3 
DO  15  J  =  1,3 

TS(I,J)  =  S*T(I,J) 
15      CONTINUE 

20  CONTINUE 

*  FIND  THE  EXPONENTIAL  OF  MATRIX  TS 

* 

CALL  EXPON(TS,EXPTS,PI) 
A  =  0.0 
DO  21  I  =  1,3 
DO  22  J  =  1,3 

A  =  A  +  AL(I)*EXPTS(I,J) 
22      CONTINUE 

21  CONTINUE 

NUMER  =  C2  -  A*(C2-C1) 
DO  30  I  =  1,3 
DO  25  J  =  1,3 

P(I,J)  =  EXPTS(I,J)  -  IDEN(I,J) 
25      CONTINUE 
30  CONTINUE 
DO  45  I  =  1,3 
DO  40  J  =  1,3 

Q(I,D  =  0.0 
D0  35K  =  1,3 

Q(I,J)  =  Q(I,J)  +Tia,K)*P(K,J) 
35  CONTINUE 

40      CONTINUE 
45  CONTINUE 
DENOM  =  0.0 
DO  55  I  =  1,3 
DO  50  J  =  1,3 

DENOM  =  DENOM  +  AL(I)*Q(I,J) 
50      CONTINUE 
55  CONTINUE 

CS(N)  =  NUMER/DENOM 
S  =  S+0.01 
IF(S.LE.L)GOTO  100 

*  FIND  THE  TIME  THAT  HAS  MINIMUM  COST 

* 
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K  =  1 
200  IF(K.EQ.N)THEN 

SMALL  =  K*0.01 
ELSEIF(CS(K+ 1),GT.CS(K))THEN 

SMALL  =  K*0.01 
ELSFJF(CS(K+ 1).LE.CS(K))THEN 

SMALL  =  (K+l)*0.01 

K  =  K  +  1 

GOTO  200 

ENDIF 

* 

*  OUTPUT  OPTIMAL  REPLACEMENT  AGE 

OPEN(UNIT  =  3,FILE  =  'OPTAGE') 

WRITE(3,333)SMALL  ,  PI 
333  FORMAT(lX,F6.3,4X,E11.4) 

OPEN(UNIT  =  4,FILE  =  'OPTCOST') 

WRITE(4,*)CS(K) 
555  STOP 

END 

SUBROUTINE  EXPON(A,EXPA,PI) 
C 

C    TfflS  SUBROUTINE  USES  TO  FIND  THE  EXPONENTIAL  OF  A  3X3  MATRIX. 
C    IT  WORK  WITH  COUPLE  OF  IMSL  SUBROUTINES 
C     ('EVCRG7LINCG7EPIRG'). 
C 

INTEGER  LDA,LDEVEC,N,NOUT,I,J,K 

PARAMETER(N=3,LDA=N,LDEVEC=N,LDUINV=N,LD=N,LDU=N) 

REAL  A(LDA,N),PI,EXPA(3,3) 

COMPLEX  EVAL(N),  EVEC(LDEVEC,N),  ELD(3,3),  UINV(3,3), 

COMPLEX  C(3,3),  D(3,3) 

OPEN(UNIT=l,FILE= 'EXPONENT') 

*  CALCULATE  THE  EIGENVALUES  AND  FJGENVECTERS 

CALL  EVCRG(N,A,LDA,EVAL,EVEC,LDEVEC) 
DO  15  I  =  1,3 
DO  10  J  =  1,3 

ELD(I,J)  =  0.0 
10      CONTINUE 
15  CONTINUE 
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ELD(U)  =  CEXP(EVAL(1)) 
ELD(2,2)  =  CEXP(EVAL(2)) 
ELD(3,3)  =  CEXP(EVAL(3)) 

*  CALCULATE  THE  INVERSE  EKjENVECTERS'  MATRIX 

CALL  LINCG(N,EVEC,LDU,UINV,LDUINV) 
DO  30  I  =  1,3 
DO  25  J  =  1,3 
C(I,J)  =  0.0 
DO20K  =  1,3 

C(I,J)  =  C(I,J)  +  ELD(I,K)*UINV(K,J) 
20  CONTINUE 

25      CONTINUE 
30  CONTINUE 
DO  45  I  =  1,3 
DO  40  J  =  1,3 
D(I,J)  =  0.0 
DO  35  K  =  1,3 

D(I,J)  =  D(I,J)  +  EVEC(I,K)*C(K,J) 
35  CONTINUE 

40      CONTINUE 
45  CONTINUE 
DO  55  I  =  1,3 
DO  50  J  =  1,3 

EXPA(I,J)  =  REAL(D(I,J)) 
50      CONTINUE 
55  CONTINUE 
C 

C  CALCULATE  THE  PERFORMANCE  INDEX  OF  EIGENVALUES  AND 
C  EIGENVECTERS  IF  LESS  THAN  1  'EXCELLENT, IF  IT  IS  BETWEEN  1 
C  AND  100  'GOOD',  IF  IT  IS  GREATER  THAN  100  'THE  SOLUTION  IS 
C  SUSPECTED' 
C 

PI  =  EPIRG(N,N,A,LDA,EVAL,EVEC,LDEVEC) 

RETURN 

END 
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PROGRAM  SIMLIF 
C 

C  THIS  PROGRAM  PROVIDES  A  RANDOM  NUMBER  OF  A  PHASE  TYPE 

C  DISTRIBUTION  THAT  HAVE  SPECIFIC  REPRESENTATION  WITH  AN 

C  INITIAL  PROBABILITIES  VECTOR(ALPHA)  AND  A  MATRIX  OF 

C  TRANSITION  RATES  BETWEEN  TRANSIENT  STATES.  THIS  PROGRAMS 

C  WORKS  WITH  SUBROUTINE  'RANNUM'  AND  MORE  SPECIFIC  CASE 

C  ONLY  4  STATES  DISTRIBUTION.  THE  VARIABLES  ARE  AS  FOLLOWS: 
C  LD(I, J)  =  TRANSITION  RATE  FROM  STATE  I  TO  STATE  J 

C  P(I,J)    =  TRANSITION  PROBABILITY  FROM  STATE  I  TO  STATE  J 

C  AL(I)    =  INITIAL  PROBABILITY  OF  BEING  IN  STATE  I 

C  ELIF     =  THE  SUPPOSED  UPPER  BOUND  EXPECTED  LIFETIME 

* 

*  INITIALIZATION 

* 

INTEGER  S,SEED,N,N1 

REAL  U,EXP,UF,ELIF 

REALLD12,LD13,LD14,LD21,LD23,LD24,LD31,LD32,LD34 

REALP12,P13,P21,P23,P31,P32,AL1,AL2,AL3,D 

DATALD12,LD13,LD14,LD21,LD23,LD24,LD31,LD32,LD34 

+      /5.0,5.0,5.0,6.0,6.0,6.0,7.0,7.0,7.0/ 

DATAP12,P13,P21,P23,P31,P32,AL1,AL2,AL3,ELIF 

+       /0.8,0.2,0.167,0.833,0.143,0.286,1.0,0.0,0.0,10.0/ 

CALL  EXCMSCFILEDEF  10  DISK  FILE  LFTEST  (DISP  MOD') 

OPEN(UNIT  =  1,FILE  ='LIFETIME') 

OPEN(UNIT  =  2,FILE  ='SEED') 

READ(2,*)SEED 

CLOSE(2) 

LIF  =  0.0 

* 

*  GENERATE  UNIFORM  0  AND  1  FOR  INITIAL  PROBABILITY 

CALL  RANNUM(1,SEED,0.0,1.0,0,U) 
IF  (U.LE.AL1)  THEN 

S  =  1 
ELSEIF(U.LE.  AL1  +  AL2)THEN 

S  =  2 
ELSE 

S  =  3 

ENDIF 

* 

*  SIMULATE  LIFETIME  WHEN  IS  IN  STATE  1 

* 
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5  CONTINUE 

IF  ((S.EQ.l)  .AND.  (LIF.LT.ELIF))  THEN 
CALL  RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P12+P13)THEN 

S  =  4 

CALL  RANNUM(3,SEED,LD14,0,0,EXP) 

LIF  =  LEF  +  EXP 
ELSEIF(U.  GT.P12)THEN 

S  =  3 

CALL  RANNUM(3,SEED,LD13,0,0,EXP) 

LIF  =  LIF  +  EXP 
ELSE 

S  =  2 

CALL  RANNUM(3,SEED,LD12,0,0,EXP) 

LIF  =  LIF  +  EXP 
ENDIF 
ENDIF 

*  SIMULATE  LIFETIME  WHEN  IS  IN  STATE  2 

IF((S.EQ.2)  .AND.  (LIF.LT.ELIF))THEN 
CALL  RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P21 +P23)THEN 

S  =  4 

CALL  RANNUM(3,SEED,LD24,0,0,EXP) 

LIF  =  LIF  +  EXP 
ELSEIF(U.GT.P21)THEN 

S  =  3 

CALL  RANNUM(3,SEED,LD23,0,0,EXP) 

LIF  =  LIF  +  EXP 
ELSE 

S  =  1 

CALL  RANNUM(3,SEED,LD21,0,0,EXP) 

LIF  =  LIF  +  EXP 
ENDIF 
ENDIF 

*  SIMULATION  LIFETIME  WHEN  IS  IN  STATE  3 

* 

IF((S.EQ.3)  .AND.  (LIF.LT.ELIF))THEN 
CALL  RANNUM(1,SEED,0.0,1.0,0,U) 
IF(U.GT.P31  +P32)THEN 
S  =  4 
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CALL  RANNUM(3,SEED,LD34,0,0,EXP) 

LIF  =  LIF  +  EXP 
ELSEIF(U.GT.P31)THEN 

S  =2 

CALL  RANNUM(3,SEED,LD32,0,0,EXP) 

LIF  =  LIF  +  EXP 
ELSE 

S  =  1 

CALL  RANNUM(3,SEED,LD31,0,0,EXP) 

LIF  =  LEF  +  EXP 
ENDIF 
ENDIF 

IF((S.LT.4)  .AND.  (LIF.LT.ELIF))GOTO  5 

* 

*  OUTPUT  TO  'FILE  LIFETIME  A' 

WRITE(10,*)LIF 

WRITE(1,*)LIF 

OPEN(UNIT  =  2,FBLE  ='SEED') 

WRITE(2,*)SEED 

STOP 

END 

SUBROUTINE  RANNUM(DISTN,SEED,RPARM1,RPARM2,IPARM,X) 
C 

C  TfflS  SUBROUTINE  PROVIDES  A  RANDOM  NUMBER  OF  7  DISTRIBUTION. 

C  IT  INTERFACES  WITH  THE  LLRANDOMH  ROUTINES  PROVIDED  IN  THE 

C  NONIMSL  LIBRARY.  THE  PARAMETER  AND  CALLING  PROCEDURE  ARE 

C  AS  FOLLOWS: 

C  DIST  =  DISTRIBUTION  TYPE  YOU  WANT  TO  SELECT  AN  INTEGER 
C  BETWEEN  1  AND  7. 

C  SEED  =  THE  RANDOM  NUMBER  SEED  YOU  WISH  TO  USE. 

C  RPARM1,RPARM2,  ANDIPARM  ARE  REAL  AND  INTEGER  PARAMETERS. 

C  PASSED  TO  THE  ROUTINE  WITH  MEANINGS  WHICH  VARY  WITH  THE 

C  TYPE  OF  DISTRIBUTION  YOU  DESIRE. 

C  X  =  THE  RETURNED  RANDOM  NUMBER,  IT  IS  ALWAYS  REAL. 

C  DISTRIBUTION  NUMBERS  AND  THE  ASSOCIATED  PARM  DEFDSHTIONS 

C  1  =  UNIFORM  ON  THE  INTERVAL  RPARM1  TO  RPARM2 

C  2  =  NORMAL  WITH  MEAN  RPARM1  AND  VARIANCE  RPARM2 

C  3  =  EXPONENTIAL  WITH  RATE  RPARM1 

C  4  =  COUCHY  WITH  A  =  RPARM1  AND  B  =  RPARM2 

C  5  =  GAMMA  WITH  SHAPE  RPARM2  AND  RATE  RPARM1 
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C      6  =  POISSON  WITH  RATE  RPARM1 
C      7  =  GEOMETRIC  WITH  P  =  RPARM1 
C 

REAL  RPARM1,RPARM2,X 

INTEGER  DISTN,SEED,IPARM,N 

REAL  TEMP,VARIAT(1) 

IF(DISTN.LE.O  .OR.  DISTN.GE.8)THEN 
WRITE(10,*)'DISTN  DOES  NOT  EXIT' 
STOP 

ENDIF 


* 


GOTO  (10,20,30,40,50,60,70),DISTN 


*  GENERATE  A  UNIFORM  BETWEEN  RPARM1  AND  RPARM2 

* 

10  CONTINUE 

IF(RPARM1-RPARM2.EQ.0)THEN 

WRITE(10,*)'ILLEGAL  EQUAL  RPARMS  IN  RANNUM' 

STOP 
ENDIF 
IF(RPARM1  .GT.RPARM2)THEN 

TEMP  =  RPARM1 

RPARM1  =  RPARM2 

RPARM2  =  TEMP 
ENDIF 

CALL  LRND(SEED, VARIAT,  1,1,0) 
VARIAT(l)  =  RPARM1  +  (RPARM2-RPARM1)*VARIAT(1) 

GOTO  99 

* 

*  GENERATE  A  NORMAL  WITH  MEAN  RPARM1  AND  STD.  DEV.  RPARM2 

20  CALL  LNORM(SEED, VARIAT,  1,1,0) 
WRITE(10,*)'NORMAL(0,1)   ',VARIAT(1) 
VARIAT(l)  =  (VARIAT(1)*RPARM2)  +  RPARM1 

GOTO  99 

* 

*  GENERATE  AN  EXPONENTIAL  WITH  RATE  RPARM1  (1/MEAN) 

* 

30  CONTINUE 

IF(RPARM1  .EQ.0.0)THEN 

WRITE(10,*)'ILLEGAL  ZERO  RATE  IN  RANNUM' 
STOP 
ENDIF 
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CALL  LEXPN(SEED,VARIAT,  1,1,0) 
VARIAT(l)  =  VARIAT(1)/RPARM1 

GOTO  99 

* 

*  GENERATE  A  COUCHY  WITH  A  =  RPARM1  AND  B  =  RPARM2 

40  CONTINUE 

IF(RPARM2.LE.0.0)THEN 

WRTTE(10,*)'ILLEGAL  COUCHY  SPREAD  IN  RANNUM,  B=  \RPARM2 

STOP 
ENDIF 

CALL  LCCHY(SEED,VARIAT,  1,1,0) 
VARIAT(l)  =  (VARIAT(1)*RPARM2)  +  RPARM1 

GOTO  99 

* 

*  GENERATE  A  GAMMA  WITH  SHAPE  RPARM2  AND  RATE  RPARM1 

50  CONTINUE 

IF(RPARM1  .LE.0.0)THEN 

WRITE(10,*)'ILLEGAL  NONPOSITWE  GAMMA  RATE  IN  RANNUM' 

STOP 
ENDIF 
IF(RPARM2.LE.0.0)THEN 

WRrTE(10,*)'ILLEGAL  SHAPE  PARAMETER  IN  RANNUM' 

STOP 
ENDIF 

CALL  LGAMA(SEED,VARIAT,1,1,0,RPARM2) 
VARLAT(l)  =  VARIAT(1)*(1.0/RPARM1) 

GOTO  99 

* 

*  GENERATE  A  POISSON  WITH  RATE  RPARM1 

* 

60  CONTINUE 

IF(RPARM1  .LE.0.0)THEN 

WRITE(10,*)'ILLEGAL  POISSON  RATE  IN  RANNUM' 

STOP 
ENDIF 
CALL  LPOIS(SEED,VARIAT,1,1,0,RPARM1) 

GOTO  99 

* 

*  GENERATE  A  GEOMETRIC  WITH  P  =  RPARM1 

* 

70  CONTINUE 
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IF(RPARM1  .LE.O.O)THEN 

WRITE(10,*)'ILLEGAL  GEOM.  PROB.  IN  RANNUM' 

STOP 
ENDIF 

CALL  LGEOM(SEED,VARIAT,1,1,0,RPARM1) 
GOTO  99 

99  CONTINUE 
X  =  VARIAT(l) 
RETURN 
END 
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PROGRAM  MOMENT 
C 

C       TfflS  PROGRAM  USES  TO  COMPUTE  THE  SECOND  AND  THIRD 
C       STANDARDIZED  MOMENTS  THAT  SUPPORTS  MOMENT  MATCHING  IN 
C       GAMS  PROGRAM. 

*  INITIALIZE  DATA 

INTEGER  LIMn\N,Nl,I,J,COUNT 

REAL  UF,  NEW(2),  REPAGE,  A,  B,  C,  EX,  EXSQ,  EXC,  EX3 

REAL  STD,  SECOND,  THIRD,  CI,  C2,  AAC,  OBJ 

REAL  RAWDAT(1000,2),  FB(1000),  P(1000),  TOT,  TCOST 

DATAC1,C2,TOT,TCOST,AAC/100.0,500.0,0.0,0.0,0.0/ 

OPEN(UNIT  =  1,FILE  ='LIFEnME\  STATUS  ='OLD') 

OPEN(UNTT  =  2,FILE  ='OPTAGE\  STATUS  ='OLD') 

OPEN(UNIT  =  3,FILE  ='RAWDATA', STATUS  ='OLD') 

OPEN(UNIT  =  4,FILE  = 'MOMENT') 

OPEN(UNIT  =  5,F1LE  ='  COUNT',  STATUS  ='OLD') 

READ(1,*)UF 

READ(2,*)REPAGE 

READ(5,*)COUNT 

PRINT*,  'COUNT= ',  COUNT 

Nl  =  25  +  COUNT 

DO  5  I  =  1,1000 

P(I)  =  0.0 

FB(I)  =  0.0 
5  CONTINUE 
DO  44  I  =  1,N1-1 

READ(3, 1 1)RAWDAT(1, 1),RAWDAT(I,2) 
1 1       F0RMAT(1X,F7.4,4X,F3. 1) 
44  CONTINUE 
CLOSE(3) 

OPEN(UNIT=3,FILE= 'RAWDATA') 

* 

*  CALCULATION 

* 

IF(LIF.LE.REPAGE)THEN 

NEW(l)  =  UF 

NEW(2)  =  1.0 
ELSE 

NEW(l)  =  REPAGE 

NEW(2)  =  0.0 
ENDIF 
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CALL  INSERT(N1,NEW,RAWDAT) 
DO  55  I  =  1,N1 

WRITE(3,33)  RAWDAT(I,1),RAWDAT(I,2) 
33  F0RMAT(1X,F7.4,4X,F3.1) 

55  CONTINUE 

IF(RAWDAT(1,2)  .NE.  0.0)THEN 

FB(1)  =  (N1-1.0)/N1 
ELSE 

FB(1)  =  1.0 
ENDIF 
DO  10  I  =  2,N1 

IF(RAWDAT(I,2).EQ.0.0)THEN 

FB(I)  =  FB(I-1) 
ELSEIFa.LT.Nl)THEN 

FBa)  =  ((N1-I)/(N1  +  1.0-I))*FB(I-1) 
ELSE 

FB(N1)  =  0.0 
ENDIF 
10  CONTINUE 

P(l)  =  l.O-FB(l) 
DO  22  I  =  2,N1 

P(I)  =  FB(I-1)-FB(I) 
22  CONTINUE 
EX  =  0.0 
EXSQ  =  0.0 
EXC    =  0.0 
EX3  =  0.0 
DO  15  J  =  1,N1 

TOT  =  TOT  +  RAWDAT(J,1) 

TCOST  =  TCOST  +  C2*RAWDAT(J,2)  +  C1*(1-RAWDAT(J,2)) 

A  =  RAWDAT(J,1)*P(J) 

B  =  (RAWDAT(J,1)**2)*P(J) 

C  =  (RAWDAT(J,1)**3)*P(J) 

EX  =  EX  +  A 

EXSQ  =  EXSQ  +  B 

EXC  =  EXC  +  C 
15  CONTINUE 

AAC  =  TCOST/TOT 

STD  =  (EXSQ  -  EX**2)**0.5 

EX3  =  EXC  +  (2*EX**3)  -  (3*EX*EXSQ) 

SECOND  =  STD/EX 

THIRD  =  EX3/(STD**3) 
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He 

*  OUTPUT 

WRTTE(4, 100)SECOND,THIRD,  AAC,EX 
100FORMAT(2X,F11.8,4X,F11.8,4X,F8.3,4X,F6.3) 

*  CLOSE(5) 

*  OPEN(UNTT  =  5,FILE  =' COUNT',  STATUS  =  'OLD') 

*  WRITE(5,*)COUNT+l 
STOP 

END 

******************************************************************** 

SUBROUTINE  INSERT(LIMIT,NEW,X) 
C       TfflS  SUBROUTINE  USES  TO  INSERT  A  PAIR  'NEW'  DATA  INTO  THE 
C       DATA  FILE  'X'  IN  ASCENDING  ORDER. 
INTEGER  UMIT,J,K 
REAL  NEW(2),X(1000,2) 
LOGICAL  DONE 
DONE  =  .FALSE. 
J  =  1 
5  IF((J.LT.LIMIT).AND.(.NOT.DONE))THEN 
IF(X(J,1).LT.NEW(1))THEN 

J  =J+1 
ELSE 

DONE  =  .TRUE. 
ENDIF 
GOTO  5 
ENDIF 

IF(J.EQ.LIMIT)THEN 
X(J,1)  =  NEW(l) 
X(J,2)  =  NEW(2) 
ELSE 

IF(J.LT.LIMIT)THEN 

DO  10  K  =  LIMIT, J+ 1,-1 
X(K,1)  =  X(K-1,1) 
X(K,2)  =  X(K-1,2) 
10  CONTINUE 

X(J,1)  =  NEW(l) 
X(J,2)  =  NEW(2) 
ENDIF 
ENDIF 
RETURN 
END 
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PROGRAM  DATA 
C 

C  THIS  PROGRAM  PROVIDES  THE  RESULT  SEQUENTIAL  ESTIMATED 

C  ENTRIES  TRANSITION  RATES  T(I,J)  AND  INITIAL  STATE  PROB. 

C  ALPAH(I),AND  ALSO  ADJUST,  PREPARE  ALL  THE  INITIAL  DATA  THAT 

C  WILL  BE  USED  FOR  THE  NEXT  SEQUENTIAL  ESTIMATION. 

C  THE  VARIABLE  ARE  AS  FOLLOWS: 

C  A,B,C 

C  D,E,F     =  THE  ESTIMATED  TRANSITION  MATPJX 

C  G,H,I 

C  AL1,AL2,AL3  =  THE  ESTIMATED  INITIAL  STATES  PROBABILITIES 

C  X    =  RANDOM  VARIABLE  OF  LIFETIME 

C  OPTAGE  =  OPTIMAL  AGE  REPLACEMENT 

C  SECOND  =  THE  SECOND  STANDARDIZED  MOMENT 

C  THIRD  =  THE  THIRD  STANDARDIZED  MOMENT 

C         OBJ  =  OBJECTIVE  FUNCTION  VALUE  OF  GAMS  PROGRAM 

* 

*      INPUT 

* 

INTEGER  COUNT 

REAL  A,B,C,D,E,F,G,H,I,AL1,AL2,AL3,X 
REALEX,OPTCS,OPTAGE,Z,SECOND,THIRD,DI,OBJ,AAC,PI 
OPEN(UNIT  =  10,FTLE  = 'OBJ', STATUS  ='OLD') 
OPEN(UNTT  =  11, FILE  =' ALPHA',  STATUS  ='OLD') 
OPEN(UNIT  =  12,FTLE  =  'BRAVO',  STATUS  ='OLD') 
OPEN(UNIT  =  13,FTLE  ='CHAREE',  STATUS  ='OLD') 
OPEN(UNIT  =  14,FTLE  = 'DELTA', STATUS  ='OLD') 
OPEN(UNIT  =  15,FTLE  ='  ECHO',  STATUS  ='OLD') 
OPEN(UNIT  =  16,FTLE  =' FOXTROT',  STATUS  ='OLD') 
OPEN(UNTT  =  17,FTLE  =' GOLF',  STATUS  ='OLD') 
OPEN(UNIT  =  18,FTLE  = 'HOTEL', STATUS  ='OLD') 
OPEN(UNIT  =  19,FTLE  = 'INDIA', STATUS  ='OLD') 
OPEN(UNIT  =  20,FTLE  ='AL1', STATUS  ='OLD') 
OPEN(UNIT  =  21,FTLE  ='AL2', STATUS  ='OLD') 
OPEN(UNIT  =  22,FTLE  ='AL3',STATUS  ='OLD') 
OPEN(UNIT  =  23,FTLE  =  ' OPTAGE',  STATUS  ='OLD') 
OPEN(UNIT  =  24,FTLE  = 'LIFETIME', STATUS  ='OLD') 
OPEN(UNIT  =  25,FTLE  = 'MOMENT',  STATUS  ='OLD') 
OPEN(UNTT  =  26,FTLE  ='  COUNT',  STATUS  ='OLD') 
OPEN(UNIT  =  27,FTLE  ='OPTCOST',  STATUS  ='OLD') 
CALL  EXCMS('FTLEDEF  50  DISK  FILE  ESTDAT1  (DISP  MOD') 
CALL  EXCMSCFTLEDEF  60  DISK  FILE  ESTDAT2  (DISP  MOD') 
CALL  EXCMS('FILEDEF  70  DISK  FILE  ESTDAT3  (DISP  MOD') 
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OPEN(UNIT  =  30,FILE  ='SCALAR2') 

READ(10,*)OBJ 

IF(OBJ  .GT.  0.05)THEN 

A  =  -1.0 

B  =  0.3 

C  =  0.3 

D  =  0.3 

E  =  -1.0 

F  =  0.3 

G  =  0.3 

H  =  0.3 

I  =  -1.0 
ELSE 

READ(11,*)A 

READ(12,*)B 

READ(13,*)C 

READ(14,*)D 

READ(15,*)E 

READ(16,*)F 

READ(17,*)G 

READ(18,*)H 

READ(19,*)I 
ENDIF 

READ(20,*)AL1 
READ(21,*)AL2 
READ(22,*)AL3 
READ(23 ,  *)OPTAGE,PI 
READ(24,*)X 

READ(25,*)SECOND,THIRD,AAC,EX 
READ(26,*)COUNT 
READ(27,*)OPTCS 
CLOSE(26) 
PPJNT*,'OBJ=',OBJ 

IF(OPTAGE  .GT.  X)THEN 

Z  =  X 

DI  =  1.0 
ELSE 

Z  =  OPTAGE 

DI  =  0.0 
ENDIF 
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*  OUTPUT 


OPEN(UNIT  =  26,FILE  =' COUNT',  STATUS  ='OLD') 

WRITE(26,*)COUNT+l 

WRTTE(50,200)COUNT,A,B,C,D,E,F,G,H,I 

WRITE(60,300)COUNT,OPTCS,OPTAGE,X,Z,DI,AAC 

WRITE(70,100)COUNT,AL1,AL2,AL3,SECOND,THIRD,OBJ,PI 


*  PREPARE  THE  DATA  FOR  "GAMS"  PROGRAM 

* 

WRITE(30,*)'     SCALAR' 

WRrTE(30,400)'  AA  INITIAL  SOLUTION  OF  A  /',A,'  /' 

WRITE(30,400)'  BB  INITIAL  SOLUTION  OF  B  /',B,'  /' 

WRTTE(30,400)'  CC  INITIAL  SOLUTION  OF  C /',C,' /' 

WRITE(30,400)'  DD  INITIAL  SOLUTION  OF  D  /',D,'  /' 

WRTTE(30,400)'  EE  INITIAL  SOLUTION  OF  A  /',E,'  /' 

WRTTE(30,400)'  FF  ESflTIAL  SOLUTION  OF  B /',F,' /' 

WRITE(30,400)'  GG  INITIAL  SOLUTION  OF  C  /',G,'  /' 

WRTTE(30,400)'  HH  INITIAL  SOLUTION  OF  D  /',H,'  /' 

WRTTE(30,400)'  H  INITIAL  SOLUTION  OF  D  /',I,'  /' 

WRITE(30,500)'  ALP1   INITIAL  SOLUTION  OF  AL1  /',AL1,'  /' 

WRTTE(30,500)'  ALP2  INITIAL  SOLUTION  OF  AL2 /',AL2,' /' 

WRITE(30,500)'  ALP3   INITIAL  SOLUTION  OF  AL3  /',AL3,'  /' 

WRITE(30,600)'  MM2  SECOND  STANDARD  MOMENT  /',SECOND,7' 

WRITE(30,600)'  MM3  THIRD   STANDARD  MOMENT /', THIRD,'/' 

WRITE(30,700)'  EX  EXPECTED  VALUE  OF  LIFETIME  /',EX,'  /  ;' 

100FORMAT(1X,I3,3(2X,F5.2),2(2X,F8.5),2X,F4.2,2X,E11.4) 

200  F0RMAT(1X,I3,2X,9(2X,F6.2)) 

300FORMAT(1X,I3,2X,F8.3,2X,F5.2,2(4X,F7.4),4X,F3.1,4X,F7.3) 

400  FORMAT(A35,F6.2,A2) 

500  FORMAT(A39,F5.2,A2) 

600  F0RMAT(A36,F11.8,A1) 

700  FORMAT(A40,F6.3,A4) 
STOP 
END 
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